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The  theory  of  the  complex  Langevin  method  for  simulating  systems  de- 
fined by  complex  valued  actions  is  examined.  It  is  argued  that  for  nonsingular 
actions,  time  independence  of  the  expectation  values  is  necessary  and  sufficient 
for  proper  convergence  of  the  stochastic  process.  Two  critera  for  determining 
a posteriori  the  accuracy  of  complex  Langevin  simulations  are  given.  By  solv- 
ing a closed  set  of  equations  for  the  moments  of  a model  problem,  it  is  found 
that  the  complex  Langevin  method  is  in  principle  more  generally  valid  than 
numerical  integrations  of  the  appropriate  stochastic  differential  equations  indi- 
cate. Singular  actions  are  discussed,  and  in  particular  the  accuracy  of  standard 
numerical  algorithms  in  the  vicinity  of  the  pole  is  examined.  For  the  singular 
model  problems  studied,  ergodicity  is  seen  to  be  frequently  violated,  and  we 
find  this  to  be  a reliable  indicator  of  the  accuracy  of  such  simulations.  The 
utility  of  the  complex  Langevin  method  for  two  dimensional  lattice  fermions 
is  discussed.  A new  measure  on  the  coherent  states  defining  a path  integral 
representation  is  introduced,  and  this  measure  is  found  to  be  quite  effective 
in  resolving  numerical  difficulties  present  in  previous  Langevin  simulations. 
Finally,  some  aspects  of  stochastic  quantization  are  discussed. 


x 


CHAPTER  1 
INTRODUCTION 


One  of  the  hard  facts  of  life  as  a physicist  today  is  that  only  a trivial 
fraction  of  the  interesting  physical  models  one  would  like  to  study  admit  exact 
solutions.  Worse  still,  reliable  numerical  results  often  prove  to  be  almost  as 
elusive;  indeed  experience  seems  to  indicate  a nagging  correlation  between 
the  interest  in  and  difficulty  of  solving  most  problems.  This  is  particularly 
true  for  models  which  contain  a large  number  of  coupled  degrees  of  freedom 
x = ... x Dy  For  such  models  the  quantities  one  would  like  to  compute 

are  often  expressed  in  the  form  of  multi- dimensional  integral  ‘averages’  over  a 
space  with  respect  to  some  ‘weight  function’  e~^: 

i r D 

(F)s  = ~r  / []  dxie~SF(x)  (1.1) 

1=1 


Given  our  complicated  and  highly  interacting  world,  such  expressions  are, 
of  course,  encountered  in  every  area  of  quantitative  science  that  seeks  to  de- 
scribe the  interaction  of  many  different  variables.  In  physics,  they  are  the 
fundamental  quantities  of  interest  in  statistical  mechanics,  quantum  mechan- 
ics, quantum  field  theory  and  condensed  matter  systems. 

When  the  number  of  variables  D is  large,  direct  numerical  integration  of 
(1.1)  based  on  a finite,  regular  discretization  (xi,X2,  ...x_/y)  of  the  space  vari- 
ables is  not  feasible,  since  the  number  of  functional  evaluations  grows  exponen- 
tially with  the  number  of  dimensions  oc  . To  appreciate  the  magnitude  of 
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such  an  endeavor,  consider  the  modest  problem  of  approximating  a 10  dimen- 
sional integral  as  a simple  Riemann  sum,  with  a regular  discretization  of  100 
intermediate  points  per  dimension.  Even  at  the  currently  impressive  perfor- 
mance rate  of  a few  Gigaflops  (1  Gflop  = one  billion  floating  point  operations 
per  second)  a modern  supercomputer  running  continuously  since  the  Renais- 
sance might  still  not  finish  its  task  in  our  century.  Understandably,  computing 
expression  such  as  (1.1)  is  an  extremely  labor  intensive  business,  and  the  prac- 
titioners of  such  black  arts  occupy  an  entire  subcommunity  within  physics,  as 
well  as  enormous  amounts  of  computer  resources. 

Perhaps  the  only  viable  alternative  to  a regular  discretization  is  to  sample 
phase  space  randomly  in  some  appropriate  way;  this  clever  idea  was  imple- 
mented as  early  as  1953  by  Metropolis  et  al.  [2].  Other  versions  and  algo- 
rithms have  followed,  but  the  generic  idea  behind  so-called  standard  Monte 
Carlo  (MC)  estimations  of  integrals  of  the  form  (1.1)  is  to  imagine  the  weight 

i n 

function  to  be  the  equilibrium  distribution  of  some  suitably  chosen 

stochastic  process  {x(r);r  > 0}.  The  process  is  constructed  so  that  its  proba- 
bility distribution  P(x,  r)  converges  in  stochastic  time  r to  the  function 
and  hence  expectation  values  computed  from  the  process  will  converge  to  the 
desired  integral  averages 

f P.  i r P. 

E(F(x(r)))  = J n dxiP(x,r)F(x)  jf  I J[dxi  e~S F(x)  r » 1 . 

i—  1 i—  1 

(1.2) 

When  the  process  is  also  ergodic,  as  one  always  hopes  will  be  the  case,  then 
ensemble  averages  may  be  replaced  by  the  appropriate  long  time  averages 

T 

P J drF(x(r))  ->  (F)s  , 

0 


T > 1 . 


(1-3) 
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Note  that  this  fictitious  ‘stochastic  time’  r is  independent  of  any  physical  time 
which  may  or  may  not  exist  in  the  original  problem. 

The  potential  power  of  such  methods  is  clear:  ND  functional  evaluations 
may  be  replaced  by  D stochastic  time  evolution  equations.  Roughly  speaking, 
if  M is  the  number  of  finite  time  step  iterations  required  for  the  simulated 
process  to  converge,  then  there  will  be  a drastic  savings  in  computational  effort 
provided  that  /D  M. 

Of  course,  there  are  many  possible  choices  for  such  a process.  Perhaps 
the  most  widely  used  currently  are  the  so-called  acceptance-rejection  methods, 
which  are  based  on  a discrete  jump-process  defined  by  the  transition  probability 
per  unit  time 

P(xu  X2)  = min  ^1,  e^Xl^_^X2^  . (1.4) 

That  is,  the  evolution  of  the  system  from  state  xi  to  state  X2  is  accepted  or 
rejected  based  upon  the  above  transition  probability.  Also  common  are  the 
so  called  heat  bath  and  molecular  dynamics  methods.  While  standard  Monte 
Carlo  algorithms  perform  well  for  many  problems  and  are  justifiably  popular, 

1 r* 

the  required  interpretation  of  as  a probability  is  clearly  restricted  to  real 

valued  actions  S.  This  is  indeed  a very  unfortunate  limitation,  since  physical 
problems  of  a quantum  mechanical  nature  generally  are  defined  by  a complex 
valued  action  S G C . 

To  circumvent  this  problem,  two  approaches  are  generally  taken  within 
the  context  of  these  methods.  First,  a Euclidean  space  formulation  (Wick 
rotation)  will  remove  the  factor  i from  the  quantum  mechanical  path  integral 
weight  e*/^.  While  this  is  sufficient  for  many  problems,  such  as  scalar  fields 
with  real  couplings,  many  systems  such  as  those  containing  fermions  and  gauge 
theories  may  still  be  defined  by  a complex  action.  In  such  cases  the  standard 
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approach  is  to  define  a real  probability  from  the  real  part  of  the  action  as 
e-ReS  Averages  with  respect  to  the  full  weight  e~ ^ are  then  computed  via 
the  identity 


It  then  becomes  necessary  to  compute  the  average  of  the  phase  e^m^.  While 
this  procedure  is  effective  for  many  problems,  there  are  important  examples 
for  which  it  is  known  to  fail.  In  particular,  when  the  average  phase  is  small, 
statistical  errors  become  uncontrollably  large.  This  is  unfortunately  the  case 
for  some  two-dimensional  lattice  fermion  systems,  such  as  the  Hubbard  model 
away  from  half-filling,  considered  by  many  to  be  a relevant  model  for  high 
temperature  superconductivity.  Despite  much  effort,  (see,  for  example,  ref. 
[3];  for  a more  general  discussion  see  also  ref.  [4].)  this  well  known  ‘sign 
problem’  remains  unresolved. 

Quite  apart  from  these  numerical  considerations,  it  would  be  highly  de- 
sirable to  have  a formalism  for  directly  computing  quantum  mechanical  path 
integrals  in  Minkowski  (real)  time,  something  which  MC  is  a priori  incapable 
of.  Thus,  there  are  compelling  reasons  for  investigating  alternatives  to  stan- 
dard MC  methods.  In  this  work  we  examine  a very  different  approach  to  the 
problem  of  simulating  complex  actions,  an  approach  which  uses  the  Langevin 
equation  to  simulate  the  complex  action  directly. 

The  Langevin  method  was  first  introduced  by  Parisi  and  Wu  [5]  for  quan- 
tum field  theories  as  an  alternative  quantization  procedure  to  standard  oper- 
ator and  path  integral  quantization  formalisms,  and  in  this  context  is  better 
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referred  to  as  ‘stochastic  quantization’.  In  the  Parisi-Wu  scheme,  the  appro- 
priate random  process  {x(r);  r > 0}  is  that  defined  by  a Langevin  equation 

= ~\^^dr  + dw^  • 

Although  originally  proposed  for  real  valued  actions,  Klauder  [6],  and  later 
Parisi  [7],  realized  that  since  there  is  no  formal  restriction  to  a real  valued  drift 
term,  the  Langevin  method  allows  for  a natural  generalization  to  cases  where  S 
is  complex.  Unlike  standard  MC  methods,  which  ignore  potentially  important 
physical  information  in  the  complex  part  of  the  action,  the  complex  Langevin 
method  uses  the  entire  complex  action  to  define  its  stochastic  process,  and 
hence  may  converge  directly  to  the  desired  distribution.  In  addition  to  a priori 
circumventing  any  potential  sign  problem,  such  a method  is  also  naturally  very 
appealing,  since  the  original  problem  is  inherently  complex  valued. 

Complex  Langevin  (CL)  methods  have  since  been  applied  to  a variety  of 
problems,  including  one  and  two-dimensional  fermion  models  [8],  [9],  [10], 
[11],  [12],  lattice  gauge  theories  [13],  [14],  the  quantum  Hall  effect  [15]  and 
model  actions  defined  in  Minkowski  space  [16]  with  quite  promising  results. 

Unfortunately,  from  the  beginning  it  was  clear  that  CL  would  not  be  with- 
out its  own  problems.  Since  the  Hamiltonian  describing  a complex  diffusion 
process  is  not  self-adjoint,  convergence  to  the  desired  stationary  state  is  not  a 
foregone  conclusion;  the  process  may  in  some  cases  exponentially  diverge  from 
the  stationary  state  instead  of  converging  to  it  [17],  [18],  [19]. 

However,  as  long  as  such  occurrences  of  ‘exploding  solutions’  are  encoun- 
tered relatively  rarely,  it  would  seem  to  be  a small  price  to  pay;  even  if  an 
occasional  simulation  failed  to  converge,  at  least  one  would  be  assured  of  good 
results  from  those  that  did.  Alas,  even  this  hope  seemed  to  be  dashed  as  exam- 
ples were  studied  for  which  CL  appears  to  converge  to  the  wrong  answer  [20], 
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[18],  an  absolutely  ‘deadly’  prospect.  This  was  understandably  a great  obstacle 
to  CL’s  widespread  use,  since  it  has  been  generally  neither  a -priori  feasible  nor 
a posteriori  possible  to  guarantee  good  answers.  In  spite  of  these  concerns, 
work  on  CL  has  continued,  since  it  gives  good  results  often  enough  to  make  it 
seem  worthwhile.  A promising  potential  solution  to  the  convergence  problem 
was  offered  in  references  [21]  [22],  [23]  and  [24]  by  adding  an  extra  kernel  to 
the  CL  equation.  Unfortunately,  while  quite  successful  for  simple  models,  it  is 
generally  far  from  clear  what  choice  of  kernel  is  needed  for  realistic  models. 

The  outline  of  this  thesis  is  as  follows:  In  chapter  2 we  will  exam  the  un- 
derlying theory  of  the  complex  Langevin  method  and  discuss  the  conditions 
under  which  it  may  be  valid.  As  the  main  result  of  this  chapter  we  conclude 
that  a necessary  and  sufficient  condition  for  proper  convergence  is  that  the  ex- 
pectation values  become  time  independent;  essentially,  if  the  process  converges 
at  all,  then  it  does  so  correctly.  A further  criteria  for  testing  a posteriori  the 
accuracy  of  a complex  Langevin  simulation  is  given  and  found  to  be  quite 
useful. 

These  ideas  will  be  further  elaborated  on  and  tested  on  a model  problem  in 
chapter  3.  By  solving  for  the  moments  of  this  process  exactly,  we  will  learn  that 
the  CL  method  is  surprisingly  more  generally  valid  than  numerical  solutions 
of  the  stochastic  differential  equation  would  indicate.  We  argue  that  this  is 
not  due  to  a degenerate  groundstate  of  the  complex  Langevin  process,  as  has 
been  speculated,  but  rather  arises  from  the  slow  convergence  rate  controlled 
by  a relatively  ‘weak’  eigenvalue  spectrum  of  the  relevant  differential  operator. 
We  also  introduce  the  notion  of  ‘stochastic  similarity’  for  complex  Langevin 
processes  and  its  discuss  implications. 
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In  chapter  4 we  discuss  the  particular  problem  of  complex  Langevin  for  sin- 
gular actions.  This  remains  an  important  open  problem,  and  we  are  unable  to 
prove  a general  result  analogous  to  the  nonsingular  case.  Rather,  the  accuracy 
of  standard  numerical  algorithms  in  the  vicinity  of  a pole  is  discussed.  Also 
the  suitability  of  the  Langevin  method  for  a few  sample  problems  with  singular 
actions  is  examined.  We  find  a frequent  lack  of  ergodicity  in  these  simulations, 
and  use  this  as  a reliable  measure  of  the  accuracy  of  the  simulation  result. 

The  applicability  of  the  CL  method  to  the  2d  lattice  fermions  is  explored 
in  chapter  5.  We  find  a new  measure  over  the  coherent  states  defining  a path 
integral  representation  to  solve  a stability  problem  that  had  existed  in  earlier 
simulations. 

Finally,  in  chapter  6 some  aspects  of  the  method  of  stochastic  quantization 
for  continuum  field  theories  are  discussed. 


CHAPTER  2 
GENERAL  THEORY 


2.1  Real  Langevin  Method 

We  will  begin  by  reviewing  the  Langevin  method  for  real  valued  actions. 
For  simplicity,  we  restrict  ourselves  in  the  following  discussion  to  the  one- 
dimensional case,  {x  = x},  since  the  generalization  to  higher  dimensions  will 
be  immediate.  For  problems  defined  by  a real  action  5,  the  Langevin  method 
estimates  integral  averages  of  the  form 


(F)  = jf  j dxF(x)  e~S^  (2.1.1) 

N = J dx  e~S^  (2.1.2) 

by  defining  a stochastic  process  via  a Langevin  equation  whose  unique  stable 
equilibrium  distribution  is  Thus,  ensemble  averages  (or  time  aver- 

ages, since  such  a process  is  ergodic)  computed  from  the  process  will  relax  to 
the  desired  integral  averages.  The  Langevin  equation  is  given  by 

dx(r)  = -----  - dr  + dw(r)  (2.1.3) 

2 UX\T ) 

where  we  have  defined  the  notation 


dS  _ as 
dx(r)  dx 


x=x(t) 


(2.1.4) 


and  w(r)  is  a standard  Wiener  process  with  zero  mean  and  covariance 


E(w(ti)w(t2 ))  = min(r1,T2). 


(2.1.5) 
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The  reason  this  works  is  as  follows.  As  shown  in  any  number  of  standard  texts 
on  stochastic  processes  (see  for  example  ref.  [25]  ),  the  process  (2.1.3)  is 
governed  by  a Fokker-Planck  (FP)  equation  for  the  probability  density  P(x,  r) 
dP(x,r ) 1 d f d dS(x)' 


dr 


2 dx  V dx  dx 


P(x,r) 


(2.1.6) 


where  ensemble  averages  over  the  Wiener  measure  are  equal  to  integral  averages 
over  P(x , r) 

E(F(x(r)))  = J dxP(x,T)F(x)  . (2.1.7) 

Furthermore,  a general  initial  density  P(x,0)  will  converge  in  time  to  the 
unique  stationary  state 

P(z,T)-.-ife-s<*>  . (2.1.8) 

We  may  demonstrate  this  last  statement  by  following  the  arguments  of  ref. 
[18].  Let 

G(x,r)  = P(a;,r)e55(x)  . (2.1.9) 

From  the  FP  equation  it  follows  that  G obeys  a Schrodinger  like  equation 
dG(x,r ) f 1 <92 


dr  V 2 dx 2 

where  the  potential  V is 


- F(x)  G(x,t)  = -HG 


(2.1.10) 


1 fdS\ 2 1 d2S 


n*)  = o -5T  -7 


(2.1.11) 


8 \dx  ) 4 dx 2 

The  Hamiltonian  H is  nonnegative  and  admits  a complete  set  of  eigenfunctions, 
thus  we  may  expand  solutions  G(x,t)  in  terms  of  the  eigenfunctions  of  H 


Hipn(x)  = Xnipn(x) 


(2.1.12) 


so  that 


1 i °° 

G(x,  t)  = + X!  an^n(x)e~XnT 


71=1 


(2.1.13) 
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1.  c . 

Because  e 2°  is  the  unique  zero  energy  eigenstate,  it  follows  that 

1 ic 

lim  G(x,r)  = — -re  5*  (2.1.14) 

r— >oo  Jv 

and  hence 

1 c 

lim  P(x,t)  = —,e  * . (2.1.15) 

T— KX)  jV 

The  fact  that  the  ground  state  is  unique  is  important  in  that  it  ensures  the 
very  useful  property  of  ergodicity,  namely  that  equilibrium  ensemble  averages 
may  be  estimated  by  long  time  avergages 

T 

lim  i f cItF(x(t))  = (F)  . (2.1.16) 

r-+oo  1 J 
0 

Ergodicity  also  means  independence  of  stochastic  variables  widely  separated  in 
time,  in  which  case  the  characteristic  function 

C(s)  = Efe‘f-lJ^(rA  (2.1.17) 

satisfies  the  ‘clustering  condition’ 

lim  C(s  + ru)  = C(s)C(r)  (2.1.18) 

U — lOO 

where  ru(r)  = r(r  — u).  This  property  is  also  useful,  for  example,  in  estimating 
statistical  errors  from  a single  time  series,  since  a long  time  average  may  be 
effectively  treated  as  many  independent  shorter  time  averages.  Ergodicity  will 
be  discussed  in  more  detail  in  the  following  chapters. 

2.2  Complex  Langevin 

For  complex  valued  actions  5 : R — > C,  the  complex  Langevin  prescription 
looks  superficially  similar.  When  S has  a well  defined  complex  extension  S(z), 
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as  we  will  assume  throughout  this  discussion,  one  may  introduce  the  complex 
Langevin  (CL)  equation 


i / \ IdS 
dz(T)  = — — dr  + dw{T) 


(2.2.1) 


2 dz(r) 

where  z(r)  = x(r)  + tt/(r)  G C,  but  u>(r)  G R as  before.  Traditionally  it  has 
been  assumed  that  there  is  an  associated  complex  valued  function  P : R — + C 
satisfying  a pseudo  Fokker-Planck  equation 


dP(x,r ) 1 d f d dS(x) 


dr 


2 dx  \ dx 


dx 


P{x,t) 


(2.2.2) 


which  also  reproduces  the  expectation  values  of  the  process 


E(F(z(t)))  = j dxP(x,r)F(x)  . (2.2.3) 

Then,  as  in  the  real  case,  ensemble  averages  of  the  process  may  converge  to 
the  desired  integral  averages,  since  (2.2.2)  has  the  same  stationary  solution 
oc  e~ ^(x).  This  clever  idea  contains  the  essence  of  complex  Langevin,  and 
allows  for  the  possibility  of  directly  mimicking  a complex  valued  ‘distribution’ 
e . Unfortunately,  in  practice,  there  are  several  potential  pitfalls  associated 
with  this.  The  most  obvious,  and  earliest  to  be  extensively  studied,  is  the 
problem  of  convergence  of  the  pseudo  FP  equation  (2.2.2).  Unlike  in  the  real 
case,  convergence  of  solutions  of  equation  (2.2.2)  to  the  desired  stationary  state 
is  no  longer  guaranteed  when  S is  complex,  but  rather,  as  we  have  seen  in  the 
previous  section,  depends  on  the  eigenvalue  spectrum  of  the  operator 

H = -\&  + V{x)-  (2-2-4) 


In  particular,  one  expects  that,  when  H has  a negative  real  part  of  an  eigen- 
value, solutions  to  eq.  (2.2.2)  will  diverge  exponentially  in  time  from  the 
stationary  state.  This  problem  has  been  discussed  at  great  length  in  ref.  [18]. 
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We  mention  only  in  passing  that  one  potentially  promising  solution  to  this 
problem  may  be  realized  by  noting,  as  already  Parisi  and  Wu  did  [5],  that  the 

1 C 

real  Langevin  process  having  as  a stationary  solution  is  not  unique.  In 

fact,  the  pseudo  FP  equation 


dP(x,  t ) 

dr 


( 9 dS(x)\ 
\dx  dx  ) 


P(x,t) 


(2.2.5) 


with  arbitrary  kernel  K(x),  also  has  the  desired  stationary  solution.  The  cor- 
responding Langevin  equation  is 


1 Op  O T r 

dz(r)  - --K  dr  + ■ ---  dr  + VKdw(r)  . (2.2.6) 

2 oz(t)  oz(t) 

Several  authors  [21],  [22],  [23],  [24]  have  advocated  tailoring  the  choice  of 
kernel  for  the  complex  Langevin  equation  so  that  the  spectrum  of  the  pseudo 
Fokker  Planck  Hamiltonian  gives  the  desired  convergence.  While  this  has  been 
successfully  implemented  for  simple  model  problems,  for  more  general,  realistic 
models,  it  is  unfortunately  unclear  what  choice  of  kernel  is  required. 

Fortunately,  experience  has  so  far  shown  that  such  cases  of  ‘exploding  so- 
lutions’ are  relatively  uncommon  in  simulations  of  physically  motivated  prob- 
lems. More  common,  at  least  in  some  toy  models  that  have  been  studied,  are 
the  unfortunate  but  intriguing  situations  in  which  long  time  CL  averages  ap- 
pear to  simply  converge  to  the  wrong  results  [20],  [18].  Understandably,  this 
disquieting  fact  has  been  a great  problem  for  the  acceptance  of  the  CL  method 
as  a practical  tool  for  computing  integral  averages,  since  it  has  been  generally 
neither  a priori  feasible  nor  a posteriori  possible  to  guarantee  the  accuracy  of 
a given  CL  simulation. 

In  an  attempt  to  understand  this  surprising  behavior,  one  might  begin 
by  examining  the  relevance  of  the  pseudo  FP  equation  (2.2.2).  While  the 
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pseudo  FP  description  may  have  been  inspired  by  analogy  with  the  real  valued 
process,  equation  (2.2.2)  is  of  course  a mathematical  convenience  which  derives 
its  ultimate  validity  from  the  nature  of  the  true  probability  density  P(x,y,r) 
of  the  process 

dx(r)  = Re  dr  + dw(r)  (2.2.7) 

dy(r)  = Im  (-5^)  dr  . (2.2.8) 

Indeed,  there  is  no  fundamental  principle  which  guarantees  the  existence  of  a 
complex  valued  function  P describing  a CL  process  in  twice  as  many  dimen- 
sions. In  particular,  this  requires  that  for  a given  P(x,y)  > 0,  there  exists  a 
P : R — ♦ C such  that 


■/ 


E(F(z))  = / dxdyP(x,y)F(x  + iy)  = dxP(x)F(x ) 


/■ 


(2.2.9) 


Asymptotic  Behavior 

Following  the  line  of  reasoning  in  ref.  [26]  and  [27],  we  discuss  the 
complex  Langevin  process  from  the  point  of  view  of  classical  function  theory. 
Quite  apart  from  the  concern  of  whether  a pseudo  FP  description  exists  is 
the  question  of  the  asymptotic  behavior  of  the  process  {x(r),  y(r)}  defined  by 
equations  (2.2.7)  and  (2.2.8).  This  behavior  is  determined  by  the  FP  equation 

^^l=TP{x,y,r)  (2.2.10) 

where 

A A < r\  2 

T = - Gx(x,y ) - G(x,y)—  - Hy(x,y ) - H(x,y)—  + (2.2.11) 

and  for  convenience  we  have  defined  the  functions 


G = Re 


(2.2.12) 


H = Im 
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dS 

dz(r ) 


(2.2.13) 


We  mention  that  since  5(2)  = u(x,y)  + iv(x,y ) is,  by  assumption,  analytic, 


G(x,y)  = - 


1 du(x,y ) 

2 dx  ’ 


H(x,y)  = ^ 


du(x,  y) 
dy 


(2.2.14) 


Also  note  that  the  equation  for  dy  has  a zero  diffusion  coefficient,  but  is  nev- 
ertheless a stochastic  equation  through  its  dependence  on  x . 

Let  us  first  assume  that  T has  a unique  stationary  solution  P0(x,y ) 


TP0(x,y)  = 0 


(2.2.15) 


with 

Po(x,y)>  0 (2.2.16) 

and 

J P0(x,y)dxdy  = 1 . (2.2.17) 

Then  there  is  a theorem  which  tells  us  that  any  initial  distribution  P(z,y,  0) 
will  converge  in  time  to  PQ(x,y ) [25]. 

On  the  other  hand,  if  the  zero  eigenvalue  of  T is  M- fold  degenerate  then 
there  are  M ergodic  classes.  In  this  case  we  have 

M 

lim  P(x,y,r)  = VcjPj(a;,y)  (2.2.18) 

T— ► OO  Z ' 

i=  1 

and  we  cannot  guarantee  convergence  of  the  process  to  the  desired  stationary 
state.  Unfortunately,  it  is  still  an  open  question  as  to  what  conditions  guarantee 
the  existence  of  a stationary  solution  at  all.  Note  that  the  CL  process  defined  by 
the  real  two-dimensional  equations  has  a singular  diffusion  matrix.  Therefore, 
a general  statement  on  the  existence  of  stationary  solutions  can  not  be  made 
based  on  the  nature  of  the  drift  terms. 
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Further,  one  can  not  exclude  the  possibility  that  the  solutions  TP0  — 0 
exist  only  in  the  sense  of  distributions  (weak  solutions).  If  this  is  the  case 
then  it  might  be  quite  difficult  to  find  (construct)  the  stationary  density.  It  is 
easy  to  see  that  a situation  like  this  can  occur  with  the  very  simple  example 
S(x)  — ax a € R+  , which  is  now  supposed  to  be  solved  by  complex 
Langevin.  The  Langevin  equation  reads: 

dx(r)  = —ax(r)dT  + dw(r)  (2.2.19) 

dy(r)  = — ay(r)dr  . (2.2.20) 

The  stationary  density  P0(x,y ) is  then  a weak  solution  and  can  be  formally 
given  as 

P0(x,y)  = J^e-ax26(y).  (2.2.21) 


Pseudo  Fokker-Planck  Equation  for  Polynomial  Actions 

We  return  now  to  the  question  of  the  existence  of  a suitable  pseudo  FP 
distribution  (2.2.9).  This  equation  has  a formal  solution,  for  which  there  are 
several  expressions.  Two  such  expressions  are  [22] 

P(x)  = J dy  e~iydxP(x,y)  = JdyP(x-iy,y)  (2.2.22) 


and  [19] 

P(x,  = E(eikz^)e~ikx  (2.2.23) 

77  k 

where  the  sum  in  the  last  expression  is  taken  to  be  continuous  (discrete)  if  the 
variable  x lives  in  a noncompact  (compact)  space. 

However,  often  such  expressions  are  indeed  merely  formal,  as  may  already 
be  seen  with  the  simple  example  of  a Gaussian  probability  density 

P(x,y)  = ^e-x2-ay2 

7 r 


, 0 < a . 


(2.2.24) 
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For  a > 1,  we  have 

^>  = \/^Tj^/<0_1)  <2'2-25> 

however,  for  a < 1,  P is  divergent.  Clearly,  P(x,  r)  will  be  a smooth  function 
only  if  the  probability  density  P(x,j/,r)  satisfies  rather  stringent  conditions. 
Although  we  find  it  is  difficult  to  prove  a general  result  without  several  technical 
assumptions,  here  we  discuss  an  approach  inspired  by  the  form  of  the  solution 
(2.2.23). 

Of  interest  is  the  behavior  of  the  expectation  values  E(ezkz^),  both  asymp- 
totically in  time  r and  as  \k\  — ► oo.  For  example,  from  eq.  (2.2.23)  we 

see  that  this  requires  that  the  expectation  values  E(eilcz(T ))  decrease  to  0 

‘rapidly  enough’  in  some  sense  as  \k\  — > oo.  For  polynomial  actions,  note 
that  e~S  G <S(R),  where  S is  the  Schwarz  space  of  C°°  functions  of  rapid 
decrease.  Since  the  Fourier  transform  maps  this  space  onto  itself,  a necessary 
condition  that  must  be  satisfied  for  proper  convergence  of  a CL  process  is 

lim  E(eik<Th  G <S(R)  . (2.2.26) 

r— too 

If  we  then  assume  that  there  is  a tq  for  which 

cT(k)  = E{eikzW)  e S(R)  (2.2.27) 

for  each  r > tq,  then  we  may  demonstrate  the  existence  of  a well  defined 
pseudo  FP  equation  and  its  proper  convergence  as  follows:  Applying  the  Ito 
calculus  rules,  namely 

= + ^ <2-2-28> 

and  the  independence  of  increments 


E(dw(r)F(z))  — 0 


(2.2.29) 
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for  any  nonanticipating  function  F(z),  we  get 


E 


fdF(z(r)) 


= E 


i d2F  i dF  ds 


V dr  ) \2  dz2{r)  2dz{r)  dz{r) 

Since  cT(k)  £ 5(R),  so  does  its  Fourier  transform  satisfy 


(2.2.30) 


P(x,  t)  = — [ dk  cT{k)e~ikx  £ 5(R)  (2.2.31) 

27T  J 

From  the  definition  of  P,  and  under  certain  assumptions  of  smoothness,  we 
may  also  write  the  expectation  values  of  polynomial  functions  F as 


E(F(z(t)))  = J dx  P(x,t)F(x)  . 


(2.2.32) 


Then  from  the  definition  of  P and  equation  (2.2.30)  we  may  partially  integrate 
to  obtain  the  pseudo  FP  equation  for  P. 


dP(x,  t)  _ 1 d 
dr  2 dx 


8S{x)  d_ 
dx  dx 


P{x,t)  . 


(2.2.2) 


For  completeness  we  mention  that  this  equation  has  two  independent,  sta- 
tionary solutions,  namely  the  desired  one 


Po(x)  = 5 


(2.2.33) 


as  well  as  a ‘spurious’  one 


x 

Ps(x)  = e~S^  J dy  eS^  . 


(2.2.34) 


The  spurious  solution  is  not  generally  positive  definite,  and  hence  automat- 
ically excluded  when  S is  real.  However,  for  S £ C,  one  is  already  dealing 
with  complex  quantities,  and  hence  there  is  no  a priori  reason  to  assume  such 
solutions  are  not  relevant  for  the  asymptotic  behavior  of  the  process.  Indeed, 
as  may  be  seen  by  differentiation 


Ps(x)  = 0( 


1 

dxS{x ) 


) for  |x|  — > oo 


(2.2.35) 
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and  hence  Ps{x)  is  also  normalizable.  Thus  there  is  at  first  glance  a tempting 
and  tantalizing  possibility  of  interpreting  these  spurious  solutions  as  being 
responsible  for  the  ‘improper’  convergence  of  some  CL  simulations. 

These  solutions  have,  of  course,  been  excluded  by  the  assumption  that 

^lirn^  P(x,  t)  <E  $(R)  . (2.2.36) 

However,  from  a practical  perspective  this  situation  is  not  very  satisfactory, 
since  it  is  very  difficult,  if  not  impossible,  to  numerically  determine  the  rapid 
(faster  than  any  polynomial)  decrease  of  a function. 

2.3  Correct  Convergence  of  the  Process 
Nevertheless,  we  will  now  see  that  it  is  not  strictly  necessary  for  P to  always 
exist  as  a classical  function.  As  an  illustrative  example,  consider  again  the  real 
valued  action  S(x ) = . For  this  action,  the  CL  equations  are,  in  terms  of  its 

real  and  imaginary  components 

dx{r)  = —xdr  + dw(r)  , dy{r)  = -ydt  . (2.3.1) 

The  equilibrium  distribution  of  this  process  is  given  by 

po(x,  y ) = -j=e~x2S(y)  (2.3.2) 

independent  of  the  initial  distribution.  If  the  initial  distribution  is  given  by 
(2.2.24),  with  a < 1,  then  at  a later  time  r, 

P(x,y,r)  = — e~x2e~ae2Ty2  . (2.3.3) 

Although  the  CL  method  is  completely  applicable  to  this  problem,  for  r < 
— ^lna,  P is  not  a function!  However,  P may  always  be  understood  as  a 
distribution  whose  action  on  a function  F(x)  is  defined  as 

P[P]  = J dxdyP{x,y)F{x  + iy) 


(2.3.4) 
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(Note  that  for  r = — ^lna  , P(x,t ) = S(x)).  In  such  cases  the  notation 
E(F(z))  = f dxP(x)F(x)  must  be  understood  to  be  merely  suggestive.  Like- 
wise, the  pseudo  Fokker-Planck  equation,  which  is  then  meaningful  only  as  an 
integral  equation 


Jix  F(X)^H  = jdxF(x)  + ) P(X,T) 


(2.3.5) 


must  be  taken  as  a purely  formal  expression  for  the  identity 


£,  f ^E(z(t))\  /I  d2F  1 dF  dS  \ 

V dr  ) \2^2(t)  2dz{r)  dz(r)J 

If  we  concentrate  only  on  the  stationary  solutions  of  this  equation, 

demonstrate  a much  more  general  and  practical  result  as  follows. 


(2.2.30) 
we  may 


Compact  Case 

We  first  discuss  the  case  in  which  x lives  on  a compact  space,  then  continue 
to  the  noncompact  case  x G R. 

Let  x G [a,  6],  and  the  action  have  periodicity  S(x+n\)  = S(x),  A = \b—a\. 
First,  suppose  S(x)  is  represented  by  a finite  Fourier  series 


and  define 


.Tf 


(2.3.6) 


Ak(r)  = E(elkz{-P)  . (2.3.7) 

Then  it  follows  immediately  from  (2.2.30)  that 

kiT ) . for  v ' „ , 

~ df—  = ~JAh  + — 9sgAq+k  . (2.3.8) 

<7 

By  direct  substitution,  the  unique  time  independent  solution  up  to  a normal- 
ization factor,  ^ = 0,  is  seen  to  be 

b 

A’k  = fjdx  eik*e~s  . 


(2.3.9) 
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solution  of  the  form 


To  see  that  this  solution  is  indeed  unique,  consider  a potential  alternative 
ition  of  the  form 


b 


a 


Substituting  equation  (2.3.10)  into  (2.3.8)  with  condition  = 0,  we  see  that 
Q must  satisfy  the  pseudo  Fokker-Planck  equation  as  well  as  the  boundary 
condition  Q(a)  = Q(b).  In  particular,  the  ‘spurious  solution’  does  not  satisfy 
Ps(a)  = Ps(b)  since 


Therefore,  Q(x)  a e s uniquely. 

One  may  be  concerned  as  to  what  extent  the  characterization  of  solutions 
Ak  in  the  form  (2.3.10)  is  itself  unique,  since  it  is  certainly  not  true  that  any 
infinite  sequence  of  numbers  may  be  expressed  in  this  way.  Nevertheless,  any 
finite  subsequence  Ak  for  |fc|  < 2K  may  be  written  as  an  integral  over  a smooth 
function 


b 


(2.3.11) 


a 


(2.3.12) 


k=—2K 


For  K large  enough  it  follows  that 


b 


(2.3.13) 


a 


where  T is  the  pseudo  FP  operator 


f =IJL  + JL 

2 dx  dx  dx 


2 dx  dx 


(2.3.14) 


Therefore 


b 


(2.3.15) 


a 
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Unfortunately,  we  do  not  presently  have  a convincing  argument  which  would 
bring  the  limit  to  the  right  of  the  integral  and  operator  T,  thus  proving  rigor- 
ously that  lim/^oo  QK  = jje~S . 

Now  let  S more  generally  be  represented  by  an  infinite  Fourier  series.  De- 
fine the  truncated  Fourier  series 


and 


SQ= 

\q\<Q 


Sqe 


^qx 


(2.3.16) 


Aqk(r)  = E(eife«W) 


(2.3.17) 


where  zq(t ) is  the  solution  of  the  CL  process  (2.2.1)  for  the  action  (2.3.16). 
The  unique  stationary  Asq  are  then  given  by 

b 

AQk  = JfJ  dxeikxe~SQ  . (2.3.18) 

a 

Since  e are  smooth,  bounded  functions  on  a compact  space  converging  to 
e , it  follows  from  the  dominated  convergence  theorem  that 

b 

AQk  = Ji  f 4*  eikX‘~S  ■ (2.3.19) 

a 


Noncompact  Case 

When  x G R,  we  are  interested  in  those  actions  for  which  e~s  vanishes  as 
M Typically,  S(x)  is  a polynomial.  We  may  extend  the  previous  results 

to  this  case  in  the  following  way.  The  idea  is  that  since  e~ § vanishes  rapidly  as 
|x|  — > oo,  we  may  approximate  it  arbitrarily  well  by  replacing  the  action  S(x) 
with  a smooth  periodic  function  with  periodicity  A,  where  A » 1.  That  is,  we 
define  a smooth  periodic  function  Sx{x  -f  nX)  = Sx{x)  so  that  Sx(x)  = S(x) 
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for  |x|  < A , where  A < A is  somewhat  arbitrary.  For  < |x|  < A,  S\ 
may  be  any  suitable  C°°  function  chosen  so  that  Sx(x)  is  smooth  at  \x\  = A. 

S(x) 


Fig.  2.1 

Making  the  polynomial  action  periodic 


(Note  that  if  A = A,  Sx(x)  can  not  in  general  be  differentiable  at  |a:|  = A). 
Since  Sx(x)  is  bounded,  smooth  and  periodic,  it  follows  from  the  previous 
discussion  that  if  zx(t)  is  the  solution  to  the  Langevin  process  (2.2.1)  with 
action  Sx,  and 


AAJfe(T)  = £l(e“*M) 


(2.3.20) 


then 
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A 


-A 


Since  e s is  of  rapid  decrease,  by  the  dominated  convergence  theorem  it  follows 
that 


(2.3.22) 


We  have  thus  shown  that  correct  convergence  holds  in  the  limit  of  arbi- 
trarily large  compact  spaces.  We  may  also  argue  the  same  result  for  ab  initio 
noncompact  spaces.  In  the  spirit  analogous  to  the  proceeding  discussion  one 
can  consider  the  time  evolution  of  the  expectation  values  An(r)  = E(zn(r )) 
which,  for  polynomial  actions  of  the  form 


Again,  we  may  easily  demonstrate  that  the  desired  stationary  solution  may  be 
written  as 


Also,  the  spurious  solution  may  be  exorcised  in  the  following  way.  When  x G R, 
Q{x)  must  satisfy  more  than  Q(— oo)  = Q(oo)  — 0.  By  direct  substitution  into 
eq.  (2.3.8)  it  is  seen  that  the  necessary  boundary  condition  is 


P 


(2.3.23) 


is  determined  by  the  equation 


(2.3.24) 


(2.3.25) 


(2.3.26) 


However,  as  may  seen  by  differentiation,  for  \x\  > 1, 


dxS(x ) 


(2.3.27) 


hence,  as  |x|  — > oo 
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(2.3.28) 


As  for  the  uniqueness  of  the  form  of  solution 


An  — 


J dx  xnQ(x) 


(2.3.29) 


we  may  attempt  to  argue  in  a similar  manner  as  in  the  compact  case.  Defining 
an  infinite  sequence  of  rapid  decrease  functions  {Qn(x)  G S(R)}  by 


. However,  as  in  the  compact  case,  we  cannot  show  directly  that  linqv->oc  Qn  = 
follows. 

Nevertheless,  with  these  provisions,  we  have  made  strong  arguments  to 
arrive  at  the  main  conclusion  of  this  chapter:  If  the  expectation  values  computed 
from  a CL  process  become  time  independent,  they  necessarily  converge  to  the 


n < N 


(2.3.30) 


we  will  again  find  that 


lim  / dx  TQ jg  = 0 


(2.3.31) 


correct  values. 


CHAPTER  3 
A MODEL  PROBLEM 


From  the  previous  chapter  we  have  seen  that  time  independence  of  the 


expectation  values  is  necessary  and  sufficient  for  their  proper  convergence. 
Another  useful  necessary  condition  is  given  by  the  fact  that,  by  the  Riemann- 
Lebesgue  lemma, 


lim  (e*kx\  = 0 . 

:l— >00  \ / 


|fc|-KX> 


(3.1) 


This  behavior  must  also  be  reflected  in  the  expectation  values  E(eikz^)  for 
large  r.  As  we  will  now  discuss,  in  numerical  simulations,  it  is  found  that  these 
two  criteria  are  excellent  indicators  of  the  accuracy  of  CL  simulations. 


3.1  Compact  Space  Model  Problem 
Let  us  consider  the  action 

S{0)  = -pcosd,  (3eC  ,9e[0,2ir)  . (3.1.1) 

The  CL  equation  in  this  case  is  given  by 

d9(r)  = -^  sin  9(t)cIt  + dw(r)  . (3.1.2) 

This  model  is  exactly  solvable  with 

(cos  0)  = /i  (/?)//<,(/?)  (3.1.3) 

where  Iv  is  a modified  Bessel  function  of  the  first  kind.  To  integrate  eq.  (3.1.2), 
an  explicit  Runge-Kutta  scheme  which  is  0(h 2)  accurate  in  the  time  step  h 
was  used.  For  the  general  SDE  with  additive  noise 

dx(r)  — a(x)dr  + dw(r) 
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(3.1.4) 
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the  RK  algorithm  is  [28] 

a0  = a(x0) 

ai  = a(xo  + aoh  + Vhrj)  (3.1.5) 

xh  — x0  + h(a0  + al)/2  + 'Shv  • 

In  principle,  ry  is  a normal,  random  variable  with  zero  mean  and  unit  variance: 
E(rf)  = 0 , E(ri2)  = 1.  Actually,  since  we  only  require  0(h2)  accuracy  of 
the  algorithm,  r)  need  not  necessarily  be  a Gaussian.  For  our  simulations  we 
chose  the  computationally  less  expensive  random  variable  r)  = 2\/3r  where  r is 
uniformly  distributed  in  the  interval  [— ■j,^]  [29].  To  within  statistical  error, 
this  choice  of  random  variable  gives  no  discernible  difference  to  the  Gaussian 
distribution,  and  often  saves  a considerable  amount  of  CPU  time. 

The  simulation  results  were  obtained  by  performing  a running  time  integral 
over  the  Langevin  time  r and  averaging  over  several  sample  paths: 

T 

TE (cos 9)  = ^ J dr  E(cos9(t))  T ->■  oo  (3.1.6) 

0 

where 

i P 

E(F)  = pY.Fi  (3.1.7) 

i= 1 

is  an  estimate  of  the  expectation  value  computed  from  the  simulated  process 
by  summing  over  P independent  sample  paths  i = 1 . . . P. 

Rather  than  assuming  ergodicity  for  individual  time  series,  as  is  often  done 
in  Monte  Carlo  type  studies,  statistical  errors  were  always  computed  over  the 
several  independent  sample  paths. 

The  series  of  figures  3.1  to  3.3  shows  for  two  values  of  /3  the  typical  be- 
havior of  CL  simulations  which  yield  the  correct  averages.  Figure  3.1  plots  the 
running  time  average  TE(cos9 ) vs.  the  integration  time  T,  while  figure  3.2 
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plots  only  the  expectation  value  E(cos9)  as  a function  of  time  r.  Figure  3.3 
shows  the  instantaneous  values  of  the  characteristic-like  function  E(eik for 
r 1.  These  all  exhibit  the  desired  behavior  as  expected  from  the  theoretical 
discussions  of  the  previous  chapter. 

Figures  3.4  to  3.6  illustrate  the  behavior  of  ‘incorrect’  CL  simulations  by 
plotting  the  same  averages  as  in  figures  3. 1-3.3.  We  see  that  for  these  simula- 
tions, the  ensemble  expectation  values  do  not  become  time  independent  and 
have  large  statistical  errors.  This  contrasts  with  the  impression  one  might  get 
from  viewing  merely  the  running  time  averages  in  fig.  3.4.  We  believe  this  is 
an  important  lesson  to  bear  in  mind  when  doing  Monte  Carlo  type  simulations 
in  general  ! Also  note  that  the  characteristic-like  function  E(e^z)  does  not 
fall  off  as  |fc|  — > oo.  In  all  cases,  the  distinction  between  correct  and  incorrect 
simulations  is  apparent. 


Re  TE(cose)  Re  TE(cos0) 
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Figure  3.1 

Real  part  of  the  time  average  TE  (cos  9)  vs.  time  T from  eq.  (3.1.2)  for  two  (3  values.  With 
100  independent  sample  paths  and  time  step  h — 0.005.  Solid  lines  indicate  exact  values. 


Re  E(cos0)  Re  E(cos0) 
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Figure  3.2 

Real  part  of  the  ensemble  average  E(cos0)  vs.  time  r from  eq.  (3.1.2)  for  two  0 values. 
With  100  independent  sample  paths  and  time  step  h — 0.005.  Solid  lines  indicate  exact 


values. 


Re  E(exp(ik0))  Re  E(exp(ik0)) 
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Figure  3.3 

Real  part  of  the  ensemble  average  E(elk6)  at  large  time  r from  eq.  (3.1.2)  for  two  0 values. 
With  100  independent  sample  paths  and  time  step  h = 0.005. 


(esoo)ax  ay  (esoo)ax  sh 
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Figure  3.4 

Real  part  of  the  time  average  T E (cos  9)  vs.  time  T from  eq.  (3.1.2)  for  two  (3  values.  With 
100  independent  sample  paths  and  time  step  h = 0.005.  Solid  lines  indicate  exact  values. 


Re  E(cos0)  Re  E(cos0) 
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Figure  3.5 

Real  part  of  the  ensemble  average  E(cos  6)  vs.  time  T from  eq.  (3.1.2)  for  two  0 values. 
With  100  independent  sample  paths  and  time  step  h — 0.005.  Solid  lines  indicate  exact 


values. 


Re  E(exp(ik0))  Re  E(exp(ikfl)) 
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Figure  3.6 

Real  part  of  the  ensemble  average  E(eikf))  at  large  time  r from  eq.  (3.1.2)  for  two  (3  values. 
With  100  independent  sample  paths  and  time  step  h = 0.005. 
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3.2  Comparison  With  Solution  for  Moments 

Even  while  the  preceding  numerical  results  confirm  the  theoretical  predic- 
tions of  chapter  2,  we  will  now  show  the  very  surprising  result  that,  contrary  to 
what  numerical  integration  of  the  SDE  (3.1.2)  indicates,  in  almost  all  cases  the 
CL  equation  is  indeed  appropriate  to  this  problem,  in  that  ensemble  averages 
from  the  stochastic  process  converge  to  the  desired  integral  averages. 

For  the  action  (3.1.1),  equation  (2.3.8)  becomes 

^ = ~Ak(r)  - *2*  (AM  - 4k-i)  - £ MkqAq(r)  . (3.2.1) 

<7 

This  set  of  equations  may  be  solved  numerically  in  the  following  way.  First 
note  that,  since  the  stationary  solution  is  unique,  the  behavior  of  solutions  for 
large  times  is  uniquely  determined  by  the  eigenvalue  spectrum  of  the  matrix 
Mkg,  independent  of  the  initial  values.  Hence  the  initial  values  Ak(r  = 0)  = SkQ 
may  be  chosen  (This  is  a convenient,  but  not  necessary  choice).  Furthermore, 
the  stationary  solutions  A^.  satisfy  |Aj? | — > 0 as  \k\  — > oo.  Therefore  we  may 
effectively  truncate  the  infinite  set  of  equations  (3.2.1)  to  a finite  set  —K  < 
k < K and  evolve  these  forward  in  time.  Of  course,  in  practical  calculations 
it  is  necessary  to  take  K large  enough  so  that  the  final  results  Ak(r  -4  large) 
become  independent  of  the  truncation.  Note  that  this  procedure  is  not  unlike 
the  method  of  coupled  Green’s  functions. 

Alternatively,  we  may  simply  solve  for  the  eigenvalues  and  zero  eigenvec- 
tor of  the  (truncated)  matrix  Mkg.  If  there  are  no  positive  eigenvalues,  the 
stationary  solution  (zero  eigenvector)  is  stable,  otherwise,  solutions  diverge. 
This  was  confirmed  with  the  EISPACK  routine  CG  [30],  and  is  completely 
consistent  with  the  numerical  solutions  to  (3.2.1). 
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Table  1 compares  the  results  of  solving  equation  (3.2.1)  with  exact  results 
(cos0)  and  results  from  the  numerical  solutions  of  the  SDE  (3.1.2).  Surpris- 
ingly, we  see  that  CL  averages  converge  to  the  proper  values  more  generally 
than  integration  of  (3.1.2)  suggests. 
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Table  1:  Comparison  of  (cos#) 

0 

arg{0 ) 

(cos0)  (exact) 

DETERMINISTIC 

SDE 

1.0 

0° 

(0.4464,0.0000) 

(0.4464,  0.0000) 

(0.4374,0.0000) 

±(0.0095,0.0000) 

30° 

(0.4255,0.1933) 

(0.4255,0.1933) 

(0.5181,0.1215) 

±(0.0107,0.0008) 

60° 

(0.3165,0.4225) 

(0.3164,0.4225) 

(0.4234,0.0709) 

±(0.0153,0.0077) 

90° 

(0.0000,0.5751) 

(0.0000,0.5749) 

no  convergence 

3.0 

0° 

(0.8100,0.0000) 

(0.8100,0.0000) 

(0.8063,0.0000) 

±(0.0032,0.0000) 

30° 

(0.8501,0.1110) 

(0.8501,0.1110) 

(0.8345,0.0972) 

±(0.0023,0.0020) 

60° 

(1.0200,0.1125) 

(1.0204,0.1125) 

(0.8369,0.1052) 

±(0.0057,0.0030) 

90° 

(0.0000,-1.3038) 

no  convergence 

no  convergence 

5.0 

0° 

(0.8934,0.0000) 

(0.8934,0.0000) 

(0.8895,0.0000) 

±(0.0012,0.0000) 

30° 

(0.9116,0.0555) 

(0.9116,0.0555) 

(0.9084,0.0557) 

±(0.0010,0.0096) 

45° 

(0.9292,0.0751) 

(0.9292,0.0751) 

(0.9245,0.0774) 

±(0.0015,0.0014) 

60° 

(0.9440,0.0996) 

no  convergence 

no  convergence  | 

90° 

(0.0000,1.8445) 

no  convergence 

noconvergence 
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To  see  if  this  behavior  is  simply  an  artifact  of  a poorly  chosen  numerical  al- 
gorithm, we  have  also  examined  the  following  numerical  schemes  for  integrating 
eq.  (3.1.2): 

Trapezoidal  Rule 

This  is  an  implicit  equation  which  must  be  iteratively  solved  at  each  time 
step.  While  computationally  more  expensive  than  the  explicit  RK  scheme,  in 
some  cases  it  is  more  stable. 

a0  = a(x0) 

«l  = a(xh)  (3.2.2) 

xh  = ^0  + Ha0  + ai)/2  + Vhr]  . 

Dynamical  Time  Step 

In  the  limit  of  small  time  step  h,  the  diffusion  term  ( O(vTi)  ) must  dom- 
inate the  drift.  In  numerical  solutions,  artificially  large  drift  terms  may  be 
avoided  by  choosing  h at  each  iteration  as 

h = min  (hmax,  1 /a2(x0))  (3.2.3) 

Reparameterizations 

One  may  also  solve  alternative  SDE’s  arising  from  a reparameterization  of 
the  variable  6 to  Cartesian  coordinates  x , y.  This  may  be  done  at  the  level  of 
the  integral  averages  or  of  the  stochastic  equations  themselves.  First,  we  may 
rewrite  the  measure  of  the  integral  averages  by  multiplying  and  dividing  by  an 
additional  factor 

00 

J drr  e~G ^ 

0 


(3.2.4) 


38 


where  G(r)  is  some  suitably  chosen  function  (e.g.  a polynomial).  Transforming 
from  polar  coordinates  (r,  0)  to  Cartesian  coordinates  ( x , y)  defines  the  new 
action 

S(x,y)  = G(x2  + y2)  - fdx/r  . (3.2.5) 

For  the  simple  choice  of  G = -^^(p+l),  this  leads  to  the  SDE’s 

dx(T)  = (fr  ~ “ Xr2P)  dT  + dWl  (T)  (3-2.6) 

dy (r)  = ~ yr2P)  dT  + dw2{r)  ■ (3.2.7) 

Alternatively,  we  may  transform  the  equation  (2.2.1)  directly  with  the  def- 
initions 

x — cosd,  y = sin#  (3.2.8) 

which,  through  the  Ito  rules  gives  the  SDE’s  with  multiplicative  noise 

dx(T ) = ^{y/x2  + 0y2  - x)dr  - ydw(r)  (3.2.9) 

dy(r)  = -(y  + ^/3xy)dr  + xdw(r)  . (3.2.10) 

After  extensive  study  we  find  that  all  of  the  above  algorithms  yielded  es- 
sentially the  same  results.  Clearly,  the  tenacity  with  which  these  simulations 
fail  requires  consideration. 


3.3  A Degenerate  Groundstate? 

While  we  have  demonstrated  that  the  ensemble  expectation  values  from 
the  process  (3.1.2)  will  generally  converge  correctly,  recall  from  the  previous 
chapter  that  the  groundstate  of  a complex  Langevin  process  may  not  be  unique. 
Hence  the  assumption  of  ergodicity  may  become  questionable  in  some  cases. 
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In  particular,  the  assumption  of  ergodicity  is  invoked  in  order  to  justify 
estimation  of  ensemble  averages  with  long  time  averages 

T 

T^oc  Y / drF^(r))  = E(F'>  ■ (3.3.1) 

0 

We  stress  that  this  issue  is  independent  of  the  potential  existence  of  a 
degenerate  (spurious)  stationary  solution  of  the  complex  valued  pseudo  FP 
equation.  In  general,  ergodicity  may  be  lost  if  there  are  two  or  more  stationary 
solutions  P{(x,  y)  which  act  to  separate  the  space  (x,y)  explored  by  sample 
paths  into  disjoint  regions.  How  might  this  arise  from  the  perspective  of  the 
complex  Langevin  equation?  Considering  the  diffusion  process  (x(r)},  we  see 
that  given  a smooth  drift  term,  the  real  valued  Brownian  motion  w(t)  should 
preclude  the  existence  of  any  ‘forbidden’  regions  of  zero  probability  along  the 
x axis.  For  compact  spaces,  we  may  make  this  statement  quite  concrete  by  the 
following  theorem  [31]:  Let  a and  b represent  the  drift  and  diffusion  coefficients 
of  a (real)  d dimensional,  homogeneous  diffusion  process,  respectively.  Both  a 
and  b are  assumed  smooth  and  b is  bounded.  If  there  exists  a 7 such  that 

x-a(x)  < 7|x|2  (3.3.2) 

holds  for  all  x in  the  complement  of  some  compact  subspace,  then  the  process 
is  ergodic.  Essentially  this  means  that  any  ‘reasonable’  processes  on  a compact 
space  will  be  ergodic. 

For  the  process  {y(r)},  however,  the  situation  is  potentially  different.  Note 
that  the  complex  Langevin  equation  (2.2.1)  really  describes  a one-dimensional 
diffusion  process  {x(r)},  since  the  stochastic  variable  y(r)  is  completely  deter- 
mined in  terms  of  x{r).  Depending  on  the  range  of  the  solution  y[r)  = y(x(r)) 
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for  the  domain  x{t ) G Ft,  it  is  not  difficult  to  imagine  a lack  of  ergodicity  along 
the  y axis. 

Furthermore,  the  stochastic  process  which  may  be  described  by  a pseudo 
FP  density  P(x)  is  not  unique.  This  may  be  seen,  for  example,  by  noting  that 
the  formal  solution  of  eq.  (2.2.9)  may  be  written  as 

P(z)  = ■—  j dk  e~lkxQ(k,  - ik ) (3.3.3) 

where 

Q(k,q)  = J dxdy  P{x,y)elkx+iqy  . (3.3.4) 

Clearly  the  probability  density  P(x , y)  is  not  uniquely  determined  by  Q(k,  —ik), 
hence  we  see  that  far  from  having  a unique  solution,  equation  (2.2.9)  has  in- 
finitely many  ! This  is  not  surprising,  since  P(x)  only  describes  the  ensemble 
averages  of  functions  of  the  special  variable  2 = x + iy.  This  fact  will  be- 
come clear  from  the  existence  of  ‘stochastically  similar  solutions’,  which  will 
be  introduced  shortly.) 

Examining  the  long  time  averages  for  each  of  many  different  sample  paths, 
we  see  that  for  those  /3  which  give  the  correct  convergence,  these  various  av- 
erages converge  nicely  to  the  same  values.  However,  for  simulations  yielding 
incorrect  results,  different  sample  paths  had  widely  varying  long  time  averages. 
This  is  seen  in  figures  3.7  and  3.8  which  are  scatter  plots  of  the  real  part  versus 
the  imaginary  part  of  T(cos2  = 6)  for  many  independent  sample  paths.  It  is 
clear,  then,  that  such  processes  are  not  ergodic  on  the  time  scales  studied. 


T(cos0)  Im  T(cos0) 
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•6  .8  1 1.2 

Re  T(cos0) 


Figure  3.7 

Scatter  plot  of  the  average  T (cos  0)  at  large  time  T from  eq.  (3.1.2)  for  two  (3  values.  With 
100  independent  sample  paths  and  time  step  h = 0.005. 


(0soo)i  III!  (esoo)i 
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•2  .4  .6  .8 

Re  T(cos0) 


Figure  3.8 


Scatter  plot  of  the  average  T(cos  6)  at  large  time  T from  eq.  (3.1.2)  for  two  (3  values.  With 
100  independent  sample  paths  and  time  step  h = 0.005. 
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Stochastic  Similarity 

However,  if  a degenerate  groundstate  is  the  source  of  this  difficulty,  there 
is  a simple  remedy  available.  To  see  this,  let  us  first  introduce  the  notion  of 
stochastic  similarity. 

We  define  a process  z(t ) G C to  be  stochastically  similar  to  a process 
z(r)  G C if  z(r)  is  not  necessarily  stochastically  equivalent  to  the  process  z(r), 
but  nevertheless  satisfies 

E(F(z(r)))  = E(F(z(t)))  (3.3.5) 

for  analytic  functions  F = F(x  + iy). 

For  any  process  z(r)  satisfying 

dz(r)  = a(z,  t)cIt  + dw(r)  (3.3.6) 

we  may  construct  a similar  process  z(r ) via 

dz(r)  = a(z,  r)dr  + dw(r)  (3.3.7) 

by  defining  a complex  valued  Wiener  process  w as 

w(t)  = Aiwi(r)  + A2w2(t)  (3.3.8) 

where  w\,W2  are  independent,  real  Wiener  processes,  and  Ai,A2  G C satisfy 
A?  + A2  = 1. 

Such  a process  w(t)  satisfies  the  usual  Ito  calculus  rules  dw2  — dr  and 
E(F dw)  = 0 for  nonanticipating  F.  These  rules  are  necessary  and  sufficient 
for  deriving  the  set  of  equations  (2.2.30),  and  hence  we  conclude  that,  if  2(0)  = 
2(0), 


E(F(z(t)))  = E(F(z(t)))  . 


(3.3.9) 
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We  learn  that  instead  of  the  complex  Langevin  process  (2.2.1),  we  may 
study  a similar  process  having  diffusion  in  both  the  x and  y directions.  The  re- 
sulting nonsingular  diffusion  matrix  should  be  sufficient  to  guarantee  ergodicity 
in  the  complex  plane.  However,  by  studying  similar  processes  for  the  /3  cos  9 
action,  for  many  values  of  X\ , we  find  essentially  the  same  results  as  before.  We 
have  also  examined  histograms  of  the  frequency  distribution  for  the  real  and 
imaginary  parts  of  z{r ) = 9(t)  for  several  different  time  values.  While  such  an 
analysis  is  inexact,  we  find  no  evidence  of  a degenerate  stationary  distribution. 
This  may  be  seen  in  figures  3.9  through  3.12,  which  show  the  histograms  at 
four  different  times  for  the  values  \/3\  = 1,  arg  = 30°  and  \(3\  = 5,  arg  = 30°. 
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Figure  3.9 

Histogram  of  Re(z)  at  different  times  T for  \(5\  = 1.0,  arg  — 30°.  1000  sample  paths 
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Im(z) 


Im(z) 


Im(z)  Im(z) 


Figure  3.10 

Histogram  of  Im(z ) at  different  times  r for  \f3\  — 1.0,  arg  = 30°.  1000  sample  paths 
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Re(z)  Re(z) 


Re(z)  Re(z) 


Figure  3.11 


Histogram  of  Re(z)  at  different  times  T for  \/3\  = 5.0,  arg  = 30°.  1000  sample  paths 
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Figure  3.12 

Histogram  of  Im(z)  at  different  times  T for  \(3\  — 5.0,  arg  = 30°.  1000  sample  paths 
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Examining  the  eigenvalue  spectrum  of  the  deterministic  equations,  it  is 
found  that  the  minimal  eigenvalue,  which  ideally  determines  the  rate  of  conver- 
gence, varies  strongly  with  \(3\.  These  values  are  plotted  versus  \(3\  for  arg  = 0° 
and  arg  = 60°  in  figure  3.13.  In  particular,  we  see  that  the  minimal  eigenvalue 
for  (3  = 1 to  be  ~ —0.5,  while  for  (3  = 3.0  or  greater,  the  eigenvalue  is  ss  —0.9  or 
greater.  We  thus  believe  that  the  source  of  the  difficulty  in  the  convergence  lies 
in  an  unfortunate  eigenvalue  spectrum,  which,  although  in  principle  adequate, 
nevertheless  causes  numerical  difficulties.  That  is,  when  the  absolute  value  of 
the  minimal  eigenvalue  is  too  small,  the  process  converges  very  slowly.  This 
conclusion  is  greatly  supported  by  the  fact  that  even  for  purely  real  (3  ~ 1,  the 
process  gives  an  approximately  correct  answer  only  after  a rather  long  time. 
Unfortunately,  in  numerically  integrating  the  Langevin  equation  for  complex 
/3,  the  unbounded  drift  oc  will  often  eventually  lead  to  an  unbounded  sam- 
ple path.  Thus,  the  simulation  may  become  unstable  before  it  can  converge 
properly. 

It  is  interesting  to  note  that  for  the  value  \/3\  = 5.0,  arg((3 ) = 60°,  the 
CL  simulation  can  appear  convergent  for  reasonable  lengths  of  time  T before 
suddenly  diverging.  Indeed,  during  these  times  the  simulations  yield  an  accu- 
rate estimate  of  the  correct  value  (cos#).  Examining  the  eigenvalues  of  M ^ 
for  this  case,  we  see  that  there  is  one  small  positive  real  part  of  an  eigenvalue 
— 0.47.  With  a fortuitous  choice  of  initial  distribution,  this  eigenvalue  may  be 
suppressed  for  even  moderately  large  times,  and  the  expectation  values  will  be 
approximately  stationary.  From  the  perspective  of  numerical  simulations,  this 
will  be  sufficient  to  yield  accurate  estimates  of  the  expectation  values. 


real  part  of  minimal  eigenvalue  real  part  of  minimal  eigenvalue 
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Figure  3.13 

Minimal  eigenvalue  versus  \ f3\ 


CHAPTER  4 
SINGULAR  ACTIONS 


When  the  distribution  P(x)  = e~^(x)  to  be  simulated  vanishes  for  some 
values  of  x,  the  corresponding  drift  term  in  the  Langevin  equation  will  be 
singular  at  these  points.  Unfortunately,  there  is  presently  no  general  result 
on  the  validity  of  Langevin  simulations  in  such  cases,  and  the  question  of  the 
proper  use  of  Langevin  methods  in  these  cases  is  still  open.  The  authors  of  ref. 
[32]  have  advocated  introducing  some  form  of  analytic  continuation  of  a delta 
function  in  the  drift,  however,  it  is  not  yet  clear  how  this  may  be  implemented 
in  practice.  Our  concern  in  this  chapter  will  be  to  address  two  modest  issues: 
First,  whether  standard  numerical  integration  algorithms  for  stochastic  differ- 
ential equations  can  accurately  integrate  a (complex)  Langevin  equation  in  the 
vicinity  of  a singularity.  Secondly,  by  studying  model  problems,  whether  the 
Langevin  method  can  accurately  simulate  distributions  with  complex  valued 
singular  actions. 

At  the  same  time,  we  believe  that,  independently  of  Monte  Carlo  studies, 
the  issue  of  singular  drift  terms  in  the  Langevin  equation  is  itself  of  great 
interest  and  warrants  further  study. 


4. 1 Langevin  Equation  With  a Simple  Pole 
Reformulating  the  problem  in  terms  of  P(x),  the  Langevin  method  is  de- 
signed to  estimate  integral  averages  of  the  form 


= f dxF(x)P(x) 
f dxP(x) 


(4.1.1) 
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by  defining  a stochastic  process  via  the  Langevin  equation 

1 dP 

dx (T)  = ^pg^dr  + dw(T ) • (4.1.2) 

Of  interest  in  this  chapter  are  distributions  which  have  zeroes  in  the  com- 
plex plane,  of  the  form 

P(x)  cx(x-  a\)Pl  (x  - a2)P2...(x  - aN)PNe~Q ^ (4.1.3) 


where  Q(x)  is  a well  behaved,  possibly  complex  valued  function.  In  the  neigh- 
borhood of  one  of  the  poles,  the  appropriate  Langevin  equation  reduces  to 


dx{r ) = 


P 


2 x(t) 


dr  + dw(r) 


(4.1.4) 


where  we  have,  without  loss  of  generality,  shifted  the  pole  to  x = 0. 

Note  that,  for  real  x,  the  pole  effectively  prevents  crossing  at  x = 0,  so 
that 

signer))  = sign(ar(r  = 0)).  (4.1.5) 


This  phenomenon  is  known  in  the  literature  as  ‘segregation’,  and  effectively 
divides  the  x axis  into  disjoint  ergodic  classes. 

Although  we  do  not  know  of  a solution  to  eq.  (4.1.4)  for  a given  Wiener 
path  w(t),  for  x G R we  may  solve  this  equation  by  an  equivalent  stochastic 
process  given  by  [33] 


x(t) 


N 


p+  i 


i= 1 


(4.1.6) 


where  Wj,  i = 1 . . .p+ 1 are  any  independent,  standard  Wiener  processes  sat- 
isfying the  initial  constraint 


x(0)  = 


\ 


p+1 

x>2(°) 


i=l 


(4.1.7) 
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This  provides  us  with  an  exact  expression  for  the  distribution  P(x,r),  (see 
next  section),  and  thus,  for  real  x,  we  may  directly  test  the  accuracy  of  standard 
numerical  integration  schemes  with  exact  results.  We  found  the  forward  Euler 
numerical  integration  scheme  with  finite  time  step  h 

^1  = xo  + -^—h  4-  VTiry  (4.1.8) 

2x0 

where  r/  is  a standard  Gaussian  variable  of  zero  mean  and  unit  variance,  to  be 
quite  accurate  for  all  initial  values  xg.  To  prevent  artificially  large  drift  terms 
appearing  for  |x|  <C  1,  the  time  step  h was  chosen  at  each  iteration  to  be  small 
enough  to  ensure  that  the  diffusion  dominated  the  drift.  That  is, 

h = mm(hmax,  |p/2x0r2)  . (4.1.9) 

For  economy  of  notation,  in  what  follows  we  will  refer  to  the  maximal  time 
step  hmax  as  simply  h. 

The  accuracy  of  this  scheme  is  illustrated  in  fig.  4.1,  which  compares  the 
characteristic  function  P(cos(A;x(t)))  computed  from  the  simulation  at  time 
r = 0.3  with  the  exact  solution  from  eq.  (4.1.6). 


Re  E(cosk0) 
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k 


Fig.  4.1 

Characteristic  function  of  the  real  process  (4.1.4)  for  p — 2,  Xq  = 0.02  at  T = 0.3, 
computed  from  forward  Euler  scheme.  With  1000  independent  sample  paths  and  time  step 
h = 0.001.  Solid  line  is  exact  solution. 


More  generally,  x may  be  a complex  variable  by  virtue  of  a complex  valued 
initial  point  zq  € C,  in  which  case  we  write 


P 


2 z(t) 


dz(r)  - 


dr  + dw(r) 


(4.1.10) 
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where  z = x + iy  € C.  A general  solution  to  eq.  (4.1.10)  is  not  known. 
However,  we  may  write  approximate  solutions  in  the  regions  \x/y\  > 1 and 
\y/x\  1 as  follows: 

For  \x/y\  » 1,  we  have,  to  0(y/x ), 


V 

dx(r)  = 77-dr  + dw 
2x 


Mr)  = -§i dr 


(4.1.11) 

(4.1.12) 


Thus,  we  may  write 


p+l 


' t=i 


x{t)  = 

, r da/x2(a) 

y(r)  = voe  2J-o' 


(4.1.13) 

(4.1.14) 


For  \y/x\  > 1, 


px 

dx(r)  = — Kdr  + dw 

2 r 


d?/(r) 


(4.1.15) 

(4.1.16) 


holds  to  0(x/y).  Since  y(r)  is  a deterministic  function  of  time,  we  see  that 


x(r)  is  a time  dependent,  linear  process 


x(t)  = x0e2^o  + j dw{a)e*  L (4.1.17) 

To 

y(T)  = y/vo~  P(T  - ro)  • (4.1.18) 

This  solution  remains  valid,  of  course,  only  so  long  as  \y/x\  1.  We 

see  that  the  region  defined  by  this  condition  is  ‘unstable’  in  the  sense  that  a 
solution  initially  in  this  region  will  soon  leave  it,  unless  x(r)  = 0 exactly.  On 
the  other  hand,  the  region  defined  by  \x/y\  1 is  ‘stable’  in  the  same  sense. 

This  is  illustrated  in  fig.  4.2 
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Stable 


Unstable 


As  expected,  it  is  clear  that  in  the  region  \y/x\  3>  1,  x will  readily  change 
sign.  However,  more  generally  we  have,  writing  z = retG , 


dy{r)  = - 


(4.1.19) 


and  therefore 


/ \ ~2  r dv/r2(cr) 

y{r)  = yoe  Jt° 


(4.1.20) 
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Hence  |y|  monotonically  decreases  to  zero,  and  solutions  will  remain  in  the 
same  half  plane  for  all  times  r. 

Although  analytic  solutions  never  cross  the  x axis,  numerical  solutions  with 
a finite  time  step  h will,  of,  course,  violate  this  condition  to  some  extent.  For 
example,  in  the  Euler  scheme 

yi  = y°~^h  (4'L21) 
any  y\  will  change  sign  whenever  Xq  + ^ h . 

Thus,  instead  of  strict  segregation  of  upper  and  lower  half  planes,  numerical 
solutions  will  be  restricted  to  a ‘narrow  band’  about  the  y axis,  whose  width 
depends  on  the  time  step  h. 

A solution  valid  for  p = 1 was  found  in  ref.  [33];  in  this  case  it  may  be 
shown  that 

, cos20  , „ , 

dr  = — — dr  + cos  6dw{r)  (4.1.22) 

z r 

and  hence  we  may  identify  a new  stochastic  time 


df  = cos2  9 dr  . 


(4.1.23) 


With  the  above  definition,  r(r)  may  be  given  by 


r(r) 


(4.1.24) 


Although  this  solution  is  somewhat  implicit,  and  does  not  give  us  direct  infor- 
mation about  the  real  part  of  2,  we  learn  several  things  from  it.  First, 


£(r2(r))  = r()  + 2f 


(4.1.25) 


with  f monotonically  increasing  with  r,  and  f < r. 
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We  may  thus  make  a very  crude  estimate  of  the  rate  of  convergence  of 
1 2/|  -t  0 as 

lny(r)  oc  r . (4.1.26) 

Relations  (4.1.25)  and  (4.1.26)  are  found  to  be  accurately  reflected  in  the 
numerical  solution  to  eq.  (4.1.10)  for  this  case. 

4.2  Moment  Equations 

Instead  of  attempting  a solution  of  eq.  (4.1.10)  directly,  we  may  study  a set 
of  deterministic  (nonrandom)  equations  satisfied  by  the  moments  E(zn(T)).  In 
particular,  for  F(z,  r)  a nonanticipating  function,  we  have,  by  the  Ito  calculus 
rules, 

(4-2i) 

and  hence 

dE(zn ) 1 , 

~ dr  = 2n(n  + P ~ ) • (4-2.2) 

For  the  even  moments  this  leads  immediately  to 

E(z2(t))  = E(z2(  0))  + (p  + l)r  (4.2.3) 

E(z4(t))  = E(z4( 0))  + 2 (p  + 3)£;(22(0))r  + (p  + 3)(p  + l)r2  . (4.2.4) 

These  expressions  are  in  excellent  agreement  with  the  numerical  results,  as 
illustrated  in  fig.  4.3,  which  shows  the  simulated  moment  E(z2(t))  vs.  r. 


Re  E(zz(r)) 
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Fig.  4.3 

Real  part  of  second  moment  of  the  complex  process  (4.1.10)  for  p = 2,  computed  from  forward 
Euler  scheme.  With  1000  independent  sample  paths  and  time  step  h = 0.001.  Solid  line  is 
exact  solution. 

However,  the  odd  moments  require  more  subtle  considerations.  As  an 
illustration,  consider  the  case  p — 2.  According  to  eq.  (4.2.2),  for  p — 2,  the 
stochastic  variable  1 /z(t)  is  a martingale,  E(\/z)  = constant.  Indeed,  if  we 
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define  the  variable  u = 1 /z,  then  a straight-forward  application  of  the  Ito  rules 
implies 


However,  contrary  to  expectation,  this  result  is  completely  at  variance  with 
those  obtained  by  numerical  integration  of  the  SDE  (4.1.10)  ! This  is  clearly 
seen  in  fig.  4.4,  which  shows  E(z(r))  vs.  r for  an  initial  value  of  zo  — 0.17  + 
*0.10  Furthermore,  numerical  results  clearly  indicate  that  for  the  simulated 
process,  \/z  is  not  a martingale. 


(4.2.5) 


clearly  showing  that  u{t)  is  a martingale.  In  this  case  we  have,  according  to 
eq.  (4.2.2) 


Re  E(z(r)) 
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Fig.  4.4 

Real  part  of  first  moment  of  the  complex  process  (4.1.10)  for  p = 2,  computed  from  forward 
Euler  scheme.  With  1000  independent  sample  paths  and  time  step  h — 0.001.  Solid  line  is 
exact  solution. 

This  discrepancy  holds  for  purely  real  simulations  as  well,  but  in  these  cases 
we  know  the  numerical  solutions  to  be  quite  accurate.  This  may  be  understood 
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as  follows:  For  p = 2 we  have  the  Bessel  process 

dr(r ) = -o It  + dw(r)  . (4.2.7) 

Defining  r = a vector  of  three  independent  Wiener  processes,  r 

has  the  distribution 


P(r,r)  = (27rr)-3/2e-(r-r°)2/2r 


(4.2.8) 


and  hence  the  distribution  for  the  variable  r = yj to2  -f-  -f-  ic2  is  given  by 

P(r, r)  = ^-e-(r2+ro )/2r  sinh(rr0/T)  . (4.2.9) 

As  expected,  this  distribution  satisfies  the  Fokker-Planck  equation 


dP(r,T)  = l_d_  (d_  _2\ 

dr  2 dr  \dr  r J ' ’ ^ 


(4.2.10) 


The  answer  to  the  apparent  paradox  involving  the  time  dependence  of  odd 
moments  lies  in  the  proper  use  of  the  Ito  calculus  rules  used  in  deriving  eq. 
(4.2.2),  according  to  which 

die2  = dr  . (4.2.11) 


The  above  expression  is  of  course  merely  formal  notation  for  the  probabilistic 
statement  E(Fdw^)  = E(Fdr).  However,  note  that  while  E(rn ) exists  for 
n > —2,  E(rn)  diverges  for  higher  Laurant  moments.  Furthermore, 


(4.2.12) 


Since  from  eq.  (4.2.10),  £'(l/r3)  = oo,  the  Ito  calculus  rule  (4.2.11)  is  not 
meaningful  in  this  case.  It  is  interesting  to  speculate  on  a proper  generaliza- 
tion to  the  Ito  rules  for  such  cases,  and  such  questions  are  currently  being 
investigated.  Unfortunately,  without  a proper  generalized  rule,  we  cannot  at 
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present  determine  a set  of  closed  equations  for  the  odd  moments  for  complex 
valued  z , and  hence  cannot  directly  test  the  accuracy  of  the  numerical  integra- 
tion schemes  for  these  moments.  However,  in  light  of  all  the  available  evidence, 
we  feel  it  is  reasonable  to  deduce  that  the  numerical  integration  scheme  is  gen- 
erally accurate. 


4.3  Stochastically  Similar  Solutions 

The  deterministic  methods  of  the  previous  section  also  allow  for  a direct 
realization  of  the  process  (4.1.10),  albeit  in  the  following,  limited  sense. 

For  the  singular  process  z(t)  in  (4.1.10)  we  define  a similar  process  z(t) 


dS(T)  = 2 W)dT  + dd>^ 


(4.3.1) 


where  w = aw i + /3w 2 is  determined  so  that  z( 0)  = 2(0).  This  process  itself 
may  be  realized  by  another  similar  process  z'(t ) 


z'(r)  = 

subject  to  the  initial  constraint 


\ 


P+1 

(r) 

i—  1 


(4.3.2) 


m = *'(o)  = 


p+i 


\ 5Z  ^- (°)  • 

N *=1 


(4.3.3) 


To  see  why  z\t ) is  only  similar  to  5(r),  note  that  from  definition  (4.1.10), 
it  follows  that 


where  we  have  defined 


dz\T ) = 0 f(  ,dr  + dw\r) 
2z  (r) 


dw,  _ Ef=i  mdwj 


(4.3.4) 


V Eg 


u,? 


(4.3.5) 
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Although  w'  satisfies  ( dw ')2  — dr  as  required,  we  nevertheless  have 

( dw')*dw ' ^ ( dw)* dw  (4.3.6) 

showing  that  w'(t ) and  w(t)  are  not  equivalent  stochastic  processes.  * 


4.4  Langevin  Simulations 

As  stated  earlier,  general  theorems  on  the  applicability  of  Langevin  meth- 
ods to  singular  actions  are  lacking.  Hence,  it  is  still  necessary  to  rely  on 
experience  and  empirical  rules  for  such  cases. 

Gaussian  type  action 
For  the  distribution 

P(x)  = xPe-ax2+bx  (4.4.1) 


where  a,  b £ C,  Re(a)  > 0,  the  relevant  Langevin  equation  is 

dz{r)  = - az  + ^b'j  dr  + dw(r)  . (4.4.2) 

Note  that  we  have,  without  loss  of  generality,  set  the  pole  at  2 = 0,  since  the 
pole  may  always  be  shifted  by  a redefinition  of  b.  We  may  immediately  write 
down  a set  of  recursion  relations  for  the  moments,  at  least  for  n > 1,  as 


dE(zn) 

dr 


ln(n+p  - 1) 


E(zn~2)  - anE(zn)  + T^-E{zn~l)  . (4.4.3) 


nb 


For  the  special  case  b = 0 we  may  solve  these  equations  in  closed  form  as 

dE(z 2) 


dr 


= {p  + 1)  - 2aE(zz) 


dE(z4) 

dr 


= (2p  + 6)E(z2)  - 4aE(z4)  . 


(4.4.4) 

(4.4.5) 


The  author  thanks  John  Klauder  for  pointing  out  this  fact. 
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For  Re(a ) > 0,  the  equilibrium  values  of  these  moments  are  given  by 


E^(oc))  = <H±D 


E(z4(oo))  = 


(2p  + 6)Q+  1) 

8 a? 


(4.4.6) 

(4.4.7) 


which  correspond  to  the  desired  integral  averages.  Hence  we  have  in  principle 
established  the  validity  of  the  complex  Langevin  method  for  this  special  case. 
Although  we  cannot  solve  the  set  of  equations  (4.4.3)  for  general  b,  we  feel  it 
is  reasonable  to  expect  that  this  case  will,  in  principle,  not  be  dramatically 
different,  in  that  equilibrium  solutions  to  (4.4.3)  will  equal  the  desired  integral 
averages. 


In  actual  practice,  however,  equilibrium  expectation  values  are  estimated 
by  long  time  averages.  This  is  valid  for  ergodic  processes,  however,  segregation 
of  the  half  planes  in  the  region  \z\  <C  1 makes  this  assumption  tenuous.  Even 
for  numerical  solutions,  which  will  ‘leak’  across  the  axes  due  to  the  finite  time 
step  h , this  can  be  expected  to  be  problematic. 

In  simulations  of  equation  (4.4.2),  it  was  found  that  for  many  values  of 
a,  b £ C,  long  time  averages  of  the  first  two  moments,  averaged  over  many 
sample  paths,  converged  to  well  within  1 to  2 percent  of  the  desired  integral 
averages.  However,  only  slightly  less  common  were  cases  of  10  — 20  percent 
error.  A good,  but  rough  indicator  of  the  accuracy  of  the  simulations  was  the 
degree  to  which  the  generated  process  was  ergodic,  as  measured  by  comparing 
the  long  time  averages  for  each  independent  sample  path. 

Two  typical  cases  are  illustrated  with  figures  4.5  and  4.6,  which  are  ‘scat- 
ter plots’  of  the  long  time  averages  for  E(z)  £ C for  100  independent  sample 
paths.  Figure  4.5,  which  indicates  a relatively  high  degree  of  ergodicity,  shows 
the  results  for  one  set  of  a,  b values.  For  these  values,  the  (real  part)  of  the 


66 


simulated  second  moment  average  was  given  by  ReE(z^)  = 4.27  ± 0.02,  com- 
pared with  the  exact  value  of  Re  (2^)  = 4.30.  In  Figure  4.6,  we  see  a much 
lesser  degree  of  ergodicity,  and  in  this  case  the  simulated  average  was  found  to 
be  ReE(z^)  = 0.71  ± 0.01  compared  with  an  exact  result  of  Re  (x2)  = 0.85. 


0 


-.1 


^ ~2 
E- 

a 

~ -.3 


-.4 


-.5 


i 1 1 | 1 1 1 | — 1 — 1 — 1 — | — r 


a=0.87+i0.50 

b=2.60+il.50 


1.6  1.8  2 2.2  2.4 

Re  TE(z) 


Fig.  4.5 

Scatter  plot  of  long  time  averages  of  the  first  moment  for  the  process  5.2.  Each  dot  represents 
one  sample  path.  With  time  step  h = 0.002  and  104  iterations. 
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Fig.  4.6 

Scatter  plot  of  long  time  averages  of  the  first  moment  for  the  process  5.2.  Each  dot  represents 
one  sample  path.  With  time  step  h = 0.002  and  10^  iterations. 

Higher  polynomial  actions  We  also  studied  the  higher  order  polynomial 


action 


P(x)  = xPe-ax2+bx-lx4 


(4.4.8) 
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and  found  very  similar  behavior,  in  that  the  ‘degree  of  ergodicity’  is  a reliable 
indicator  of  the  accuracy  of  the  simulation  results.  At  the  same  time,  we 
mention  that  the  authors  of  ref.  [32]  have  concluded  that  strict  segregation 
is  not  necessarily  responsible  for  the  lack  of  egodicity,  as  has  been  previously 
thought. 


CHAPTER  5 
LATTICE  FERMIONS 

Fermions  are  notoriously  difficult  to  deal  with  in  any  numerical  context. 
Because  of  the  anticommuting  algebra  that  fermionic  operators  satisfy,  some 
creativity  must  be  applied  in  order  to  express  fermion  path  integrals  as  integral 
averages  over  c-number  variables.  One  common  trick  is  to  effectively  integrate 
out  the  fermionic  variables  from  the  path  integral  by  introducing  an  auxiliary 
bosonic  field,  known  as  a Hubbard-Stratonovich  field.  Here  we  discuss  quite 
a different  approach  developed  by  several  authors  [34],  [9],  [10],  [11],  [8], 
[12],  which  has  proven  suitable  for  complex  Langevin  simulations. 

5.1  Jordan- Wiener  Transformation 

We  begin  by  recalling  that  any  one-dimensional  ordered  chain  of  lattice 
fermion  operators  C\,  ..Cn,  satisfying  an  anticommuting  algebra 

{C^  ni  Cm}  = fimn  {CmCm}  — 0 (5.1.1) 

may  be  represented  via  the  so  called  Jordan- Wigner  transformation  [35]  given 
by 


c 1 = <rl 

(5,1.2) 

n—  1 

Cn  = en  JJ  ai  ’ 12  > 1 

i—  1 

(5.1.3) 

where  a is  the  spin  lowering  operator 

an  = + *>n)  (5.1.4) 
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and  <7* , &n ) &n  are  the  Pauli  matrices  at  lattice  site  n. 

Even  though  individual  operators  satisfy 

{<7+,cr_}  = l,  ((t_)2  = 0 (5.1.5) 

the  nonlocal  product  of  the  a ^ is  necessary  in  order  for  the  anticommuting 
relations  to  hold  at  different  lattice  sites.  With  this  transformation  any  lat- 
tice fermion  Hamiltonian  may  be  treated  as  an  equivalent  spin  Hamiltonian. 
When  the  original  lattice  fermion  system  is  one-dimensional,  local  and  nearest- 
neighbor  operator  products  take  on  a particularly  simple  form  in  the  spin  lan- 
guage: 

C^nCn  — ~(1  + &Z  n)  (5.1.6) 

CKCn+ 1 + h.c.  = <JXn<Txn+ 1 + <ry„<ryn+ 1 . (5.1.7) 

This  representation  may  be  extended  to  higher-dimensional  lattices  as 
well,  as  long  as  the  lattice  sites  are  ordered  into  a one  dimensional  chain 
n — 1 ,...,JV.  For  a two  dimensional  lattice  L x M — N,  a convenient  or- 
dering is  given  by  n — L x (m  — 1)  -f  /,  which  may  be  visualized  as  the  ‘zig-zag’ 
pattern  shown  below  on  a 4 x 4 lattice. 
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1=  1 . . .L 

Fig.  5.1 

Lattice  Ordering 

As  opposed  to  the  simply  periodic  boundary  conditions  one  might  be 
tempted  to  use 

(l,M  + 1)  = (/,  1)  (5.1.8) 

( L + 1 ,m)  = (l,m) 

we  impose,  following  reference  [8],  the  ‘twisted  torus’  boundary  conditions 

(l,M  + !)  = (/  + !,  1)  (5.1.9) 


( L + l,m)  = (1,  m) 
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so  that  the  resulting  one- dimensional  lattice  is  now  simply  periodic. 

While  operator  products  of  same-site  and  nearest-neighbor-sites  in  the  ver- 
tical (m)  direction  retain  their  simple  form  as  in  the  Id  case,  nearest-neighbor 
fermion  products  perpendicular  to  the  ordering  (/)  direction  take  on  a nonlocal 
form  in  their  Id  representation: 

^ 1-1  d"  h.c.  = nCji+i  + h.c.  — > ax nax n_j_i  + n+l  (5.1.10) 

M— 1 

C^lmCl+l,m  + h-C-  = C^nCn+M  + h.C.-*(crXnPXn+M+(TynPyn+M)  JJ  °Z  i 

i—n+ 1 

As  can  be  seen,  the  asymmetry  of  the  lattice  ordering  scheme  in  fig.  5.1 
introduces  an  asymmetry  in  the  form  of  nonlocal  operator  products.  This 
fact  was  the  source  of  the  numerical  difficulties  which  plagued  lattice  fermion 
simulations  in  more  than  one  dimension,  until  the  introduction  of  a modified 
coherent  state  representation  of  the  spin  operators,  which  we  now  discuss. 


5.2  Spin  Coherent  State  Representation 
If  we  wish  to  apply  Langevin  methods  it  is  necessary  to  have  a representa- 
tion of  the  appropriate  degrees  of  freedom  as  continuous  c-number  variables. 
This  may  be  readily  obtained  via  the  spin  1/2  coherent  states  |ft)  = \6,<f>) 
which  are  labeled  by  the  continuous  variables  6 , <j>  lying  on  the  two-sphere 
0<#<7r,  0 < <j)  < 2ir  . The  spin  1/2  coherent  states  are  defined  as 


IO\  _ (e  l<^  cos  6/2 

1 ' ~ V e^sme/2 


(5.2.1) 


These  states  form  an  overcomplete,  continuous  labeling  of  the  Hilbert  space 
and  thus,  for  a suitable  measure  provide  a resolution  of  unity 


1 = J dfi\Cl)(fl\ 


(5.2.2) 
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as  well  as  a representation  of  the  Pauli  matrices 


ax  = wx 

aV  = Wy 

oz  — wz 


J dfi\Q  )(  ft | x 

J y 


(5.2.3) 

(5.2.4) 

(5.2.5) 


where  x — sin0sin</>,  y = sinf?sin<^,  z = cos^.  The  ‘weight  factors’  wx,wy,wz 


are  given  by 


wx  = 


f dV 

J dfj,x ^ ’ 


wy  = 


I dd 
f dyy2  ’ 


W] 


fdV 

f dHZ2 


(5.2.6) 


and  thus  clearly  satisfy 

1 1 1 

— + — + — = 1 . (5.2.7) 

WX  Wy  WZ 

It  is  appropriate  to  refer  to  the  coefficients  wx,Wy,w2  as  ‘weight  factors’  since 
they  indicate  the  relative  weight  with  which  the  measure  y favors  larger  values 
of  X ,y  ,z  respectively.  It  follows  that  W{  > 1 for  any  absolutely  continuous 
measure  (i.  Indeed  for  uniform  weighting  by  wx  = wy  — wz  = 3. 

The  magic  of  coherent  states  lies  in  their  over  completeness,  or  equivalently, 
in  their  nontrivial  overlap,  given  by 


(5.2.8) 


This  implies  that  the  choice  of  measure  /i  is  not  unique,  and  indeed  my  be  cho- 
sen to  suit  the  particular  problem  at  hand.  As  we  shall  see,  it  is  precisely  this 
freedom  that  allows  us  to  ‘soften’  the  potentially  harsh  effects  of  the  nonlocal 
terms  for  the  2 d fermionic  system. 
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For  a fermionic  Hamiltonian  H , its  representation  in  terms  of  the  spin 
coherent  states  will  be  of  the  general  form 

H = Jdi,\{n))({[l)\h({Sl})  (5.2.9) 

where 

!{«}>{{«}!  = |n,  >(«,(  x ...  x |fi«)(SI„|  (5.2.10) 

and  is  some  function  of  the  variables  {17}  = . . . 0/y,  <f>]\r}.  The 

correspondence  H —*  h may  be  found  easily,  and  in  particular  we  have  for  the 
following  bilinear  products 

nCn  — y d~  wzzn ) (5.2.11) 

C^n+lCn  + h.c.  -+  ^(wlxn+ixn  + w^yn+iyn)  (5.2.12) 

1 M-i 

^ n+M^n  T h.C.  ► ^(wxxn+lxn  T ^yVn-^Wn)  1 1 ^z^n-\-i  • (5.2.13) 

i- 1 

We  see  now  that  for  a 2 d,  L x M = N lattice  with  the  ordering  (5.1.9)  , 
the  ‘nonlocal’  terms  contain  a potentially  large  factor  of  relative  to  the 

local  terms.  While  it  is  true  that  from  (5.2.6)  wz\z\  « 0(1)  ‘on  average’,  this 
is  generally  very  difficult  to  capture  numerically  if  is  much  larger  than 

1.  In  particular,  if  //  is  the  usual  spherical  measure 

d/z  = d£l  = — sin  6d9d<f>  (5.2.14) 

27T 

it  follows  from  (5.2.6)  that  wx  = wy  = wz  = 3,  and  hence  for  any  reason- 
ably sized  lattice,  the  nonlocal  term  completely  dominates  the  action,  making 
simulations  unstable  [8]. 
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From  (5.2.6)  it  is  seen  that  wz  will  be  reduced  by  a measure  that  is  weighted 
more  heavily  at  the  polar  regions  \z\  « 1.  One  particularly  useful  choice  of 
measure  was  found  in  ref.  [12]  to  be 


d[i  — f(Cl)dQ, 

= (1  - a)(l  + q)zq  + a 

where  0 < a < 1 and  q is  an  even  integer.  For  this  measure 

3(3  + 3) 


(5.2.15) 


WX  = Wy  — 


and 


wz 


(3  + aq) 
3 (q  + 3) 


(5.2.16) 


(5.2.17) 


(3  + [3  - 2a]q) 

The  free  parameters  a and  q allow  the  weight  factors  to  be  favorably  controlled. 
Of  course,  by  virtue  of  the  constraint  condition  (5.2.7)  we  may  reduce  wz  only 
at  the  expense  of  increasing  the  values  of  wx,  wy.  Note  however,  that  the  factors 
containing  wx,wy  appear  equally  in  the  local  and  nonlocal  terms.  Figure  5.2 
shows  the  values  of  the  weights  wx(=  wy ) and  wz  vs.  a for  q = 6 and  q = 10 


values. 


Weight  Factor  Weight  Factor 
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Fig.  5.2 


Weight  Factors  vs.  a 
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5.3  The  Partition  Function 

The  averages  we  are  interested  in  calculating  are  of  the  form 

Tr  (O  e~PH) 

(0)  = v J 


(5.3.1) 


Tr(e~PH) 

for  some  operator  0 and  /?  some  fixed  temperature.  Factoring  e~PH  — ( e~eH)K 
where  e = P]  K,  we  may  write,  correct  to  O(e^)  (suppressing  lattice  indices  n ) 


„ — cH 


Taking  the  trace 


= J d/i |J2)<fi|  e~€h  + 0(e2) 
Tr(-)  = J dfl  (ft|  • |0) 


(5.3.2) 


(5.3.3) 


we  have  an  approximate  expression  for  the  partition  function 
Tr(e-'J")  = J j)  dnke~s  + 0(e2K) 


(5.3.4) 


k=l 


where  the  action  is  given  by 


S =^2  -ln(^fc+ll ~ In f(£lk)  + eh(Qk)  . 
k= 1 


(5.3.5) 


Note  that  we  have  imposed  periodic  boundary  conditions  in  the  ‘thermal  time’ 
direction 

I^A'+l  )=  l^l)  • 


(5.3.6) 


The  numerator  of  (5.3.1)  takes  the  form 

I< 


Tr (Oe-fX)  = j J]  + °^K)  ' ^ 

For  purposes  of  the  simulation  we  symmetrize  the  above  expression  over  the 
thermal  time  lattice  to  obtain 


Tr  (Oe~PH)  = 
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In  applications  to  lattice  spin  and  lattice  fermion  problems  the  operator  0 
will  itself  be  expressed  in  terms  of  the  Pauli  matrices  ax,ay,az.  Although  the 
ratios  in  the  above  expressions  may  appear  complicated,  in  fact  they  reduce  to 
a simple  form  [36] 


_ sjb+1  + xk  + izk+lVk  ~ izicVk+l 
(fyfc+1  Ifyfc)  1 + xkxk+i  + ykyk+ 1 + zkzk+l 


(^ib+lkyl^t)  _ Vk+ 1 + yk  + izk+ ixk  - izkxk+l 
(ttk+l  Ifyfc)  1 + xkxk+\  + ykyk+ 1 + zkzk+ 1 


^ zk+ 1 +zk+*yk+lxk  ~iVkxk+ 1 (5  3 10) 

(^+ll^/fc)  1 + xkxk+i  + ykyk+i  + zkzk+i 

We  make  the  important  remark  that  since  the  spin  coherent  states  span 
the  two  dimensional  Hilbert  space,  the  overlap  function  sometimes 

consists  of  two  orthogonal  vectors,  and  hence  will  vanish  for  certain  values  of 
the  coordinates.  Specifically  this  happens  when  Qk,£lk+i  characterize  any  two 
antipodal  points  on  the  two-sphere.  Since  this  overlap  appears  as  a logarithm  in 
the  action  S,  this  may  be  expected  to  cause  some  difficulties  in  the  simulations, 
as  was  often  the  case  with  the  singular  actions  in  the  previous  chapter. 


5.4  Complex  Langevin  on 

When  computing  integral  averages  over  a space  with  a nonflat  metric  such 
as 

(F(0,  <!>))=  jfj  dne~SF(6,  <j>)  (5.4.1) 

some  care  must  be  taken  in  defining  the  appropriate  stochastic  process.  Since 
the  Langevin  equations  are  local  in  space,  we  define  the  process  in  terms  of  the 
two  locally  orthogonal  vectors  dO  and  sin  6d(j)  on  S2,  as  shown  in  fig.  5.3. 
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Since  the  effective  action  is  given  by  5 — In  sin#,  the  appropriate  Langevin 
equations  are 

1 d 

dO{r)  = --^—[S  - ]n(sm0)\dT  + dwi(r)  (5.4.2) 

• \ 1 d [5  — ln(sin#)]  . . 

sin (O)d(p(T)  = . - ( . — dr  + dw2(r)  (5.4.3) 

2 sin  uO(p(T) 

where  w\  ,w2  are  independent,  standard  Wiener  processes. 

For  numerical  simulations,  it  is  convenient  to  consider  the  above  problem 
in  the  embedding  space  coordinates  x,y,z.  Using  the  Ito  calculus  rules  for 
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change  of  variable,  in  ref.  [9]  it  was  found  that  equations  (5.4.2)  and  (5.4.3) 
are  stochastically  equivalent  to 


dx(r ) = Dxdr  + ydw^(r)  — zdw^r) 


(5.4.4) 


dy(r)  = Dydr  + zdw\(r ) — xdw^(r) 
dz(r)  = Dzdr  + xdw^r)  — ydw\(r) 


(5.4.5) 

(5.4.6) 


where  the  drift  terms  are  given  by 


D’  = ~x~2 


(l-z2) 


dS 


- xy 


dS 


dS 


xz- 


dx(r ) dy(r)  dz(r) 


d«  = -v--2 


Dr  = —Z  — ~ 


dS  , ..  2.  dS  dS 

~xy^77^  + (!  - y )- 


dx(r ) 


- yz 


dS  dS 

-xz--  , x — yz- 


dy{r ) dz(r). 

+ a-S)-ds 


(5.4.7) 

(5.4.8) 

dx(r)  dy(r)  ' v*  ~ J dz(r)\  (5.4.9) 

The  most  important  advantage  of  this  parameterization  is  the  removal  of 
the  unphysical  coordinate  singularities  at  8 = 0 and  9 = tt.  Also,  although 
the  above  requires  the  solution  of  three  stochastic  differential  equations  with 
multiplicative  noise  (as  opposed  to  the  original  two  equations  with  additive 
noise),  we  have  found  that  due  to  the  inordinate  amount  of  time  involved  in 
the  computation  of  trigonometric  functions,  the  Cartesian  parameterization  is 
even  less  CPU  intensive. 

Notice  that  the  diffusion  terms  are,  in  vector  notation,  given  by  the  vector 
cross  product  of  the  coordinate  vector  r = (z,y,z)  and  the  vector  Wiener  pro- 
cess d w,  ensuring  that  the  diffusion  is  tangential  to  the  spherical  surface.  This 
preserves  the  important  constraint  x2  + y2  + z2  = 1.  While  equations  (5.4.4), 
(5.4.5)  and  (5.4.6)  preserves  this  constraint  exactly,  in  numerical  integrations 
it  will  only  hold  to  finite  order,  and  hence  must  be  imposed  ‘by  hand’  after 
each  time  iteration.  For  actions  of  the  form  (5.3.5),  the  appropriate  CL  drift 
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terms  have  been  found  in  ref.  [9],  modified  with  the  new  measure  terms,  to 
be 

d* = - x"  + 55{(1  ~ - 4"1)- 

(4+1 + >4+1)44(4_1  - *4~‘)+ 

(i  - 4+1)[i  - 4(4  + ■4)](4-‘  - ;4“‘)+  (s.4.10) 

(4+1  + i4+1)[i  - 4(4  - ;4](i  - 4-1))- 
- (4)2l4  - 444  - 44<4) 

o,  = - 4 + ^{(i  - 4+1)44d  - 4-1)- 
(4+1+i4+1)44(4“1-i4‘1)+ 

(i  - 4+1)[i  - 4(4 + <4)l(4_1  - >4_1)+  (5.4.H) 

(4+1  + .•4+‘)[i  + 4<x‘  - i4i(i  - 4-1))- 
j{-444  + [i-(4)2]4-444} 

D>  = - 4 + 5j;{(i  - 4+1)44(i  - 4-1)- 
(4+1 + ;4+1)44(4-‘  - <4_1)+ 

(1  - 4+1)[l  - 4(x‘  + i'4)](x*_1  - ■4‘1)+  (5-4.12) 

(4+1  + *4+1)[i  - 4(4  - «4l(i  - 4-1)}- 
i{[i- (4)2]4-444 -44c*) 

where 

D =(x*+1 +.'4+1)(i  + 4)(x‘-1  - .'4-')+ 

(4+i4)(i-4+1)(4~1-i4~1)+ 

(5.4.13) 

(4+1  + >'4+1)(i  - 4)(4  - *4)+ 
(i-4)d-4“1)(i-4+1) 

and  the  functions  A,B,C  are  given  by  [9],  [12] 


(5.4.14) 
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r,k  _ dhk 
Bn  ~ ‘dy* 

(5.4.15) 

k dhk  1 dfi 

dzn  fn  dzn 

(5.4.16) 

As  an  aside,  we  remark  here  that,  as  an  alternative  to  the  expression  (5.3.2), 

we  may  also  take  to  0(e2) 

e~eH  = J d£l\Sl)(Sl\  e~€hf  + 0(e2) 

(5.4.17) 

which  would  lead  to  the  expressions 

<\k  _ fk  dhk 

OX* 

(5.4.18) 

T>k  _ rk  dhk 

“ tJn  o l 
^Vn 

(5.4.19) 

rk  _ rkdhk  ,k  dfn 

dzn  ozn 

(5.4.20) 

While  the  above  is  in  principle  equally  valid,  expression  (5.3.2)  has  a few 
practical  advantages  over  (5.4.17).  First,  since 


% = ,(s  + 1)",“1 

the  prefactor  in  (5.4.17)  may  become  quite  large  when  q « 10,  as  will  be 
typically  the  case.  Furthermore,  the  expression  involves  computing  the  entire 
lattice  energy  at  each  time  step,  a computationally  expensive  proposition  for 
large  lattices. 


5.5  Higher  Order  Accuracy 

We  may  improve  the  accuracy  of  the  partition  function  (5.3.4)  to  be  correct 
to  0(e2)  as  follows:  With 


h = / dfi\n){n\h(si) 


(5.5.1) 


we  may  write,  correct  to  0(e3), 
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-eH 


= J dndji\Sl )(  )( 


1 - \e(h(Sl)  + h(£l))  + le2h(tt)h(ti) 

Z z 


+ 0(e3) 

(5.5.2) 


which  may  be  exponentiated,  keeping  0(e3)  accuracy 

e~eH  = J dfidfi\n){n\\Ci)(n\  e-Hw)4c2(A4)2  + 0(e3)  (5.5.3) 


With  this  expression,  and  after  some  manipulations,  the  new  action  for  the 


partition  function  (5.3.4)  is  given  by 

K I<  K—l 

s = hk  + -e2  ^2  h2k  - e2  hkhk+l  . (5.5.4) 

k—1  fc=l  k odd 

Note  that,  as  two  coherent  states  were  used  in  the  representation  , K is  assumed 


to  be  even  in  the  above  expression.  To  this  order  in  e,  the  functions  Ak , Bk,  Ck 
in  equations  (5.4.14),  (5.4.15)  and  (5.4.16)  are  modified  to  be 


k odd 


k even 


Ak  — e^f(1  + e^k  ~ e^Jfc+l ) 

Bk  ~ eT“£(1  + e^k  ~ e^k+ 1) 
OVn 

Ck  = (!  + ehk  - thk+l) 


Ak  = + e^k  ~ e^k- 1) 

Bk  = e^f(1  + e%  ~ e^k-l) 
Ck  — e^f(1  + e^k  ~ ehk-l) 


(5.5.5) 

(5.5.6) 

(5.5.7) 

(5.5.8) 

(5.5.9) 
(5.5.10) 


This  procedure  may  in  principle  be  generalized  to  any  desired  order.  It  is  note- 


worthy that  the  use  of  coherent-state  techniques  lead  fairly  straightforwardly 
to  the  introduction  of  improved  actions  in  the  sense  of  higher  order  accuracy. 
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5.6  Numerical  Integration 

As  the  stochastic  differential  equations  (5.4.4)  through  (5.4.6)  involve  mul- 
tiplicative noise,  we  used  a numerical  integration  scheme  due  to  Petersen  [37] 
as  the  one  best  suited  for  our  lattice  fermion  simulations.  For  the  general 
multivariate  SDE  with  multiplicative  noise 

dxi  = a\x)dr  + bj(x)dw 3 (5.6.1) 

the  two-step  algorithm  is  given  by 


a0  ~ a\xo) 

(5.6.2) 

% = b){x o) 

(5.6.3) 

a\  = a*  (4  + Vhglojr)J0  + halQ) 

(5.6.4) 

b\j  = blj{x  o + y/h/2blQjriJ1  + ha^/ 2) 

(5.6.5) 

b2  j = bj(x  q — y/h/2bQjT)3^  + ha^/2) 

(5.6.6) 

b0j,k  = dkb)(x0) 

(5,6.7) 

xh  = x o + 2^ao  + ai)  + j + b2  jWo  + hbokbbi,jAkl 

(5.6.8) 

where  Tj q and  rj\  are  independent,  random,  normal  variables  with 

zero  mean 

and  unit  variance.  The  random  matrix  A lJ  is  given  by 

(5.6.9) 

and  is  of  the  form 

T”  = 1 

(5.6.10) 

Tij  = -Tji  ,i^j 

(5.6.11) 

where  each  element  T1^  is  itself  an  independent,  standard  normal  random  vari- 
able as  well.  This  algorithm  is  0(h?)  accurate  for  all  moments  E(xn). 
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From  the  above  equations,  the  diffusion  matrix  bl-  is  given  by 

/ 0 -2  y \ 

b=  2 0 -x  ) (5.6.12) 

\-y  x 0 / 

Hence,  b1-  k = — c,-^.  The  diffusion  terms  proportional  to  h in  equation  (5.6.8) 
axe  then  given  by 


bikAkl  = yA12  + zA13 


- x(A22  + ,433) 


(5.6.13) 


bh^kAkl  = xA21  + zA23  - y(An  + A33)  (5.6.14) 

4lAkAkl  = xA*1  + VA™  ~ z(aU  + A 22)  (5.6.15) 


where  we  have  suppressed  all  lattice  indices.  The  stochastic  matrix  A^  (5.6.9) 
is  given  by 


( (^o)2  ~ 1 

vfoo  + Tl2 

W^l-T13 


vWo  ~ T 12 
(^o)2  ~ 1 
vWo  + T 23 


vWo  - T 23 

(.03)2  - 1 


(5.6.16) 


5.7  A Test  Case 

We  wish  to  test  the  appropriateness  of  the  CL  method  for  computing  in- 
tegral averages  of  the  form 

i r K 

<°> = ir  / n int  ° e~s  <5-71) 

J k=l 

with  the  simple  example  of 

h(Qk)  = pzk8lk  , mk)  = 1 . (5.7.2) 

This  choice  of  h(£lk)  does  not  correspond  to  any  Hamiltonian,  but  does  provide 
an  exactly  solvable  toy  model.  In  particular,  CL  was  used  to  compute  ( zk ), 
which  is  given  by 


(zt)  = (coth/3  - 1/3}  6lk  . 


(5.7.3) 
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This  is  shown  in  fig.  5.4,  where  TE(z\)  is  plotted  against  /?. 


0 


Fig.  5.4 

TE(z\)  vs.  f3  for  K=4. 

As  can  be  seen,  there  is  excellent  agreement  with  exact  results  for  this  toy 


model. 
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5.8  Free  Fermion  Model 

We  consider  now  a two-dimensional,  nonrelativistic,  spinless  free  fermion 
model  with  a finite  chemical  potential  //,  which  is  defined  by  the  Hamiltonian 

h = -tj2c}Cj-»j2cici  <5-81) 

(u) 

where  the  sum  Yl(ij)  over  nearest  neighbors  lattice  sites. 

This  Hamiltonian  can,  of  course,  be  diagonalized  in  terms  of  momentum 
eigenstate  creation  and  annihilation  operators 

H = J2EkCkCk  (5.8.2) 

k 

where  the  energy  eigenvalues  are  given  by 


Ek  = — 2f(cos  kx  + cos  ky) 


(5.8.3) 


and  the  Bloch  wave  vectors  k = ( kx,ky ) are  restricted  to  the  first  Brillion  zone 

(5.8.4) 

(5.8.5) 


2irmx  2tt  mv 

kx  = — T7 7T  , kv  — 


— 7T 


N " ’ N 

mx,my  = 1 ...  N . 

For  free  particles  satisfying  Fermi-Dirac  statistics,  the  partition  function  in  the 
grand  canonical  ensemble  is  given  by 

Q = II  l1  + . (5.8.6) 

k 

Thus  the  average  particle  occupation  number,  or  filling,  is  given  by 


-i 


(5.8.7) 


Returning  now  to  the  complex  Langevin  algorithm,  in  this  representation, 
the  appropriate  complex  action  for  the  Hamiltonian  model  is: 


N K 

S = ^2~ln((x^y^z)n+1\^,yiZ)n)  -ln/(4) 

n=l  (fc=l 


(5.8.8) 
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-~Ywx(xn+lxn  + Vn+lVn) 

. M- 1 

_£°I  2 m-i,  k , Jfc.Jfe  TT  J 

2 wxwz  \xn+Mxn  > Vn+Myn)  zn+j 

j=  1 

a0/C,  , fcv 

2~(1  +^z2„)  • 

As  a first  study  we  are  interested  in  the  effect  that  the  choice  of  measure 
(5.2.15)  over  the  coherent  state  labels  9,<f>  has  on  the  distribution  generated 
by  the  simulations.  From  previous  studies  [8]  it  was  found  that  the  natural 
measure  on  the  sphere  dO  = ^ sin  0d6d<f)  is  insufficient  for  simulations  on  two 
dimensional  lattices  of  any  reasonable  size.  We  confirmed  this  result  and  found 
simulations  with  this  measure  to  be  highly  unstable  for  lattice  sizes  greater  than 
3x3.  Figure  5.5  shows  the  distribution  of  the  real  parts  of  the  variable  z versus 
the  real  part  of  x at  one  lattice  site  of  a 3 x 3 lattice,  summed  over  all  ensembles 
and  imaginary  time  slices  I\ . Comparing  this  with  Fig.  5.6,  which  shows  the 
same  plot  for  a 4 x 4 lattice,  it  is  clear  that  the  larger  lattice  shows  instabilities 
quickly,  in  that  the  variables  quickly  obtain  large  values  outside  the  surface  of 
the  two-sphere,  leading  to  floating  point  overflows.  The  situation  worsens  for 
larger  lattices.  Figures  5.7  and  5.8  show  the  results  of  the  simulations  with 
the  same  physical  parameters,  but  with  the  measure  (5.2.15),  and  the  choices 
of  a = 0.2,5  = 10?  and  — 0.1,5  = 8 respectively.  Clearly  there  is  marked 
improvement,  and  the  simulations  are  stabilized  quite  satisfactorily. 
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Fig.  5.5 

Distribution  of  X,  Z on  the  sphere  for  3x3  for  a = 1.0.  = 1.0  fl  — 2.0 
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Fig.  5.6 

Distribution  of  x , Z on  the  sphere  for  4 X 4 for  a = 1.0  . (3  = 1.0  pi  = 2.0 
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Fig.  5.7 

Distribution  of  x , 2 on  the  sphere  for  4x4  for  a = 0.2,  q = 10. 
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Fig.  5.8 

Distribution  of  X,  z on  the  sphere  for  4x4  for  a = 0.1,  q = 8 
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For  the  free  model  our  order  parameter  of  interest  was  the  filling 


which  as  a spin  operator  is  given  as 


(5.8.9) 


n=  1 

First  we  explored  the  sensitivity  of  the  simulation  results  to  variations  of 
the  free  parameters  a,  q.  Table  2 compares  the  values  obtained  for  T E(p)  for 
several  a values  for  q — 8 and  q = 10  respectively,  for  a 5 x 5 x 5 lattice, 
with  f3  = 1.0, /i  = 2.0.  In  principle,  of  course,  the  results  are  completely  inde- 
pendent of  the  free  parameters.  However,  as  we  have  discussed,  the  numerical 
cancellations  necessary  for  this  are  delicate  for  larger  lattices,  and  hence  we 
expect  some  variation  in  the  simulations  results  with  the  choice  of  a,  q.  As 
an  estimate  of  the  statistical  error,  we  took  the  standard  deviation  a over  all 
independent  sample  paths,  instead  of  the  conventional  statistical  error  a/\/N, 
as  the  former  reflects  more  accurately  the  variation  in  results  among  the  vari- 
ous sample  paths.  It  may  be  seen  from  Table  2 that  the  results  are  reasonably 
robust,  within  statistical  error,  with  respect  to  variations  in  a,  q. 
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Table  2:  Filling  for  various  a,  q 

<1 

a 

wz 

wx 

TE(P) 

8.0 

0.05 

1.269 

9.435 

0.82  ±0.04 

0.10 

1.299 

8.689 

0.80  ±0.07 

0.15 

1.341 

7.865 

0.86  ± 0.03 

0.20 

1.386 

7.181 

0.85  ± 0.03 

0.25 

1.435 

6.598 

0.84  ±0.04 

10.0 

0.20 

1.345 

7.797 

0.85  ±0.03 

0.25 

1.393 

7.089 

0.86  ±0.04 

0.30 

1.444 

6.504 

0.84  ±0.05 

0.35 

1.500 

6.000 

0.87  ±0.04 

0.40 

1.560 

5.571 

0.86  ±0.04 

Examining  the  time  dependence  of  the  expectation  value  E(p(r)),  we  find 
that  the  value  appears  approximately  stationary  after  3000  to  4000  iterations. 
Figure  5.9  compares  E(p)  versus  r for  the  same  parameters,  but  two  different 
initial  distributions  of  the  variables  x , y,  z.  As  can  be  seen,  there  is  a moderate, 
but  not  significant,  dependence  of  the  asymptotic  values  on  the  initial  condi- 
tions. This  may  be  seen  more  clearly  by  comparing  the  long  time  averages 
T{p)  for  each  independent  sample  path.  The  scatter  plot  in  fig.  5.10  shows  a 
noticeable  variation  in  these  averages  for  different  paths,  indicating  a lack  of 
‘strong’  ergodicity  on  the  time  scales  studied. 

We  believe  this  lack  of  ergodicity  is  due  to  the  singular  nature  of  the  co- 
herent state  overlap  functions  in  the  action,  and  is  similar  to  the  situation 
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with  singular  actions  studied  in  the  previous  chapter.  Under  the  present  cir- 
cumstances, we  find  the  general  method  described  thus  far  to  be  insufficiently 
accurate  to  reliably  obtain  exact  values  for  the  filling.  However,  it  is  also 
known  that  this  algorithm  may  give  an  accurate  qualitative  description  of  the 
behavior  of  lattice  fermions  for  various  interacting  models  [9] . We  have  found 
this  to  be  also  the  case  with  the  free  fermion  model,  as  is  illustrated  in  fig.  11, 
which  plots  the  filling  versus  the  chemical  potential  p for  a 5 x 5 x 5 lattice, 
P = 0.3. 

One  may  wonder  why  the  toy  model  of  section  5.7  seemed  to  work  so 
well.  Notice  that  for  that  choice  of  the  partition  function  factorizes  into 

independent  integrals  for  each  degree  freedom.  This  fact  is  reflected,  albeit 
indirectly,  in  the  Langevin  equations.  Given  an  accurate  numerical  scheme, 
all  degrees  of  freedom  will  decouple  in  this  way,  leading  effectively  to  a single 
degree  of  freedom  problem,  free  of  singularities. 

We  conclude  that,  while  quite  promising,  more  work  is  required  before  the 
method  we  have  described  can  be  consistently  and  accurately  applied  to  lattice 
fermion  systems. 


Re  E(p)  Re  E(p) 
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Fig.  5.9 


E(p)vs.T  for  two  different  initial  distributions.  For  10  sample  paths,  a 


0.1,g  = 8 
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Fig.  5.10 

Scatter  plot  ofT(p)  for  two  different  initial  distributions.  For  10  sample  paths,  a = 0.1,  q = 8 


98 


Fig.  5.11 

Filling  for  several  fJ,  values. 

5.9  Interacting  fermions 

While  we  have  seen  that  the  algorithm  we  have  presented  for  fermion  sim- 
ulations still  has  problems  to  be  overcome,  we  feel  that  its  great  potential  for 
studying  interacting  fermions  justifies  continued  strong  interest  in  the  method. 
One  of  the  most  promising  features  arises  from  the  fact  that  interaction  terms 
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in  the  fermion  Hamiltonian  generally  introduce  no  real  complication  to  the 
algorithm.  That  is,  free  fermions  and  interacting  fermions  appear  virtually  the 
same  from  the  standpoint  of  the  effective  action  S.  To  see  this,  consider  the 
two-dimensional  single  band  Hubbard  model,  which  is  defined  by  the  Hamil- 
tonian 

B = -*  £ ClCi,  - „ 5]  clCi„  + U £ (5.9.1) 

(ija)  i<?  * 

This  model  differs  from  the  free  Hamiltonian  (5.8.1)  by  the  introduction  of  spin 
degrees  of  freedom  a =|,  which  are  summed  over  in  the  above  expression, 
as  well  as  by  an  interaction  characterized  by  a coupling  constant  U.  Since  the 
fermion  number  operator  C^C  is  represented  (via  the  Jordan- Wigner  transfo- 
mation)  in  the  spin  operator  language  as  C^C  = ^(1  + &z),  the  contribution  of 
the  interaction  term  to  the  effective  action  in  Cartesian  coordinates  is  simply 

N K U 

hI  = +wzznO(1  +wzzh)  (5-9‘2) 

n—  1 k—l 

This  shows  that  the  introduction  of  interactions  of  this  form  in  the  Hamiltonian 
leads  only  to  the  addition  of  quadratic  terms  in  the  effective  action.  Thus,  in 
this  scheme,  interactions  are  only  a minor  modification  to  the  free  case. 


CHAPTER  6 

STOCHASTIC  QUANTIZATION 


Up  until  now  we  have  been  mainly  concerned  with  numerical  applications  of 
the  Langevin  equation.  However,  it  is  important  to  realize  that  the  present  use 
of  and  interest  in  Langevin  techniques  for  numerical  Monte  Carlo  studies  were, 
to  a large  extent,  inspired  by  their  use  as  an  alternative  quantization  scheme  for 
continuum  quantum  field  theories.  In  order  to  bring  things  more  in  perspective, 
we  now  briefly  discuss  the  beautiful  subject  of  stochastic  quantization , a field 
which  began  with  an  original  paper  by  Parisi  and  Wu  [5]  in  1981.  The  literature 
on  such  a mature  field  is  vast  and  there  are  excellent  reviews  books  on  the 
subject  [38],  [39],  [40].  In  this  chapter  we  will  simply  highlight  a few  main 
aspects  of  a very  large  subject,  which,  even  if  it  were  not  for  some  of  its  practical 
advantages  would  be  worthy  of  study  on  its  aesthetic  merits  alone. 

Of  course,  we  have  already  been  exposed  to  the  idea  of  stochastic  quanti- 
zation in  its  zero-dimensional  (finite  number  of  degrees  of  freedom)  form.  For 
continuum  quantum  fields,  the  action  of  a Euclidean  space  path  integral  is  pre- 
sumed to  be  the  equilibrium  distribution  of  an  infinite-dimensional  stochastic 
process.  The  choice  of  a Euclidean  space  formulation  is,  of  course,  designed  to 
avoid  the  complex  valued  equation  that  would  naturally  arise  from  the  factor 
i/h,  in  the  Minkowski  space  Feynman  path  integral.  Most,  but  not  all  work  in 
the  field  presumes  this  formulation.  (More  will  be  said  on  this  point  later.) 

Consider  the  simplest  case  of  a single  scalar  field  for  which 


(m)  = 


JV(j)  e~SEF{(j> ) 

f V(f>  e~^E 


(6.1) 
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where  Sg  denotes  the  Euclidean  action,  and  we  have  set  h = 1.  In  the  following 
we  will  drop  the  subscript,  and  simply  write  the  appropriate  Langevin  equation 
as 

1 6S 

d<£(x,r)  = --— dr  + dw(x,T ) (6.2) 

Z 0(p{X,  T) 

where 

E(w(x,Ti)w(y,T2))  = 8(x  - y)min(Ti,T2)  (6.3) 


and  8S/8<j>  denotes  the  functional  derivative  of  the  total  action 


S=  j d^xCi^dfi) 


(6.4) 


with  respect  to  the  field  <f>(x).  Note  that 


8S 

8<f>(x) 


(6.5) 


yields  the  classical  field  equations,  so  that,  as  expected,  the  process  spends  most 
of  the  time  varying  about  the  classical  solution.  Also,  reinserting  h we  see  that 
in  the  semi-classical  limit  h,  <C  1 the  drift  term  clearly  dominates  the  process, 
although  the  unbounded  Wiener  paths  keep  the  process  from  becoming  com- 
pletely classical.  Strictly  speaking,  since  the  variables  <f>(x)  are  distributions 
with  potentially  infinite  operator  products,  solutions  to  (6.2)  generally  exist 
only  in  the  weak  sense,  i.e.,  when  integrated  against  a suitable  test  function. 

In  deference  to  the  almost  universal  notation  and  conventions  in  use  on 
this  subject,  we  write  the  Langevin  equation  formally  as 


d<f)(x,  r) 
dr 


SS 

8(f) 


+ l{x,r) 


(6.6) 


where  r/(x,r)  are  the  independent  white  noises  satisfying 


E(ri(x,  T)7)(y,  f ))  = 26(x  - y)6(r  - f ) . 


(6.7) 
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Formally,  t](t)  = \/2^p,  however,  as  may  be  seen  by  the  delta  correlations 
(nondifferentiability  of  the  Wiener  paths),  this  process  is  not  pointwise  defined, 
that  is,  it  is  a weak  equation.  Note  that  the  factor  of  ^ missing  in  the  above  drift 
term  appears  instead  in  the  noise  correlation.  This  is  completely  equivalent  to 
the  previous  formulation,  and  has  a few  aesthetic  merits. 

Of  course,  when  discussing  continuum  field  theories  with  potentially  infinite 
renormalizations,  the  derivations  of  Chapter  1 showing  the  convergence  of  the 
expectation  values  to  the  desired  equilibrium  values 


V in/  j t \ u \ \ fv<t><j){x  i)...^(x„)e  5 

W(*i  ’ r)  • • • ’ r))  = jvfps 


(6.8) 


expectedly  take  on  a more  formal  character.  Indeed  the  concern  of  many 
earlier  studies  in  the  field  was  to  demonstrate  convincingly  that  the  stochastic 
quantization  prescription  does  indeed  yield  the  traditional  interpretation  of 
the  Feynman  path  integral.  Here  we  venture  our  personal  opinion  that,  where 
continuum  field  theories  are  concerned,  quantization  is,  to  some  extent,  in  the 
eye  of  the  beholder.  One  might  consider  with  equal  weight  the  suggestion  that 
a priori  stochastic  quantization  is  at  least  as  valid  a prescription  as  any. 

Consider  the  specific  example  of  a scalar  field  with  a traditional  looking 
action 


~Jfx 


5W)2  + + V(4.) 


(6.9) 


where  the  local  potential  may  be  typically  a polynomial  with  some  small  pa- 
rameter such  as  \<j)4. 

Just  as  in  normal  QFT,  one  is  interested  in  obtaining  the  free  form  of 
the  propagator,  and  then  performing  a perturbation  expansion  about  the  free 
solution.  * 

* There  are  compelling  objections  to  this  point  of  view  [41]. 
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The  free  theory  is  Gaussian  in  nature,  so  it  is  no  surprise  that  the  solution 
to  the  Langevin  equation  for  this  case  may  be  given  in  momentum  space  (of 
the  space-time  variables)  as 

T 

(/>{k,  t)  = ^(Jfe,  0)e-r(P+m2)  + J dt  e-tr-tHP+m2)^)  (6<10) 

0 

where 

</>(A:,r)  = J dx  e^x 4>(x,t)  (6-11) 

and  0)  is  the  initial  field  value  at  r = 0.  Hence  we  obtain  the  momentum 
space  propagator  as 

E(<Kk,  r ■)*(«,  r))  = (2*)"  (l  - (6.12) 

which  clearly  shows  that  as  r — ► oo,  we  get  the  desired  Euclidean  space  free 
propagator. 

In  a manner  quite  similar  to  standard  perturbative  QFT,  one  may  develop 
a perturbation  series  expansion  of  solutions  to  the  Langevin  equation.  That 
is,  we  consider  the  retarded  propagator  G(t)  in  k space 


dG 

dr 


— ( k '?  + m2)G  + S(r) 


(6.13) 


with  the  boundary  condition  G(k,r)  = 0 , r < 0.  The  solution  is  given  by 


OU.t)  = e-(*!+m2>T0(r)  . 


(6.14) 


Then  the  general  solution  to  the  Langevin  equation  may  be  written  as  an 
integral  equation.  For  the  example  of  V(<f>)  = \<ffi  this  is  given  by 

T 


^(jfe,r)  = J dt  c-(*2+ma)(r-0 
0 


•KM)- j%+  ^(p.<)^(g.0*(*-p-«) 


(6.15) 
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which  may  be  solved  order  by  order  in  A 

* = JG1+^J  J J (G’lM’lMn)  + . . . (6.16) 

Already  in  their  original  paper  Parisi  and  Wu  introduced  the  diagrammati- 
cal rules  similar  to  the  traditional  Feynman  rules  for  solving  the  above  equation 
order  by  order  in  A,  and  indeed,  after  an  integration  over  the  white  noise  r 7, 
an  exact  correspondence  can  be  shown  between  the  Feynman  diagrams  and 
the  ‘stochastic’  diagrams  in  the  limit  r — > 00.  We  mention  only  in  passing 
that  the  machinery  of  stochastic  quantization  is  flexible  enough  to  ‘prepare’ 
the  ensemble  at  r = —00,  leading  effectively  to  a stationary  process  for  finite 
r,  so  the  r — >•  00  limit  is  not  strictly  necessary. 

A major  impetus  behind  the  development  of  the  stochastic  quantization 
method  has  been  its  potential  for  quantizing  gauge  theories  without  the  need  to 
explicitly  fix  the  gauge,  as  is  the  case  in  the  standard  path  integral  quantization 
scheme  (Since  the  action  is  gauge  invariant,  an  integration  over  the  gauge 
degrees  of  freedom  would  lead  to  a divergent  volume  in  the  partition  function). 

For  an  Abelian  gauge  theory  defined  by  the  Maxwell  action 

s = \j  dnxF^F^  (6.17) 


the  appropriate  Langevin  equation  is  then 

dA 


dr 


= dvF^v  + ^(r) 


(6,18) 


where  the  white  noise  correlation  is  given  by 


E{rin(x,  T)rju(y , t ))  = 2 S(x  - y)S(r  - f )Sfiiy  . (6.19) 


The  Langevin  equation  is  of  course  invariant  under  gauge  transformations 
of  the  form 


Afx(x,  r)  — > Afi(x,  t ) + dfiA(x)  . 


(6.20) 
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In  momentum  space,  the  Langevin  equation  is 


dAfj 

dr 


— (fc  k^ku)Au  -f-  T)n 


(6.21) 


Gauge  invariance  of  the  action  implies  that  the  transverse  field 


— k^ki>  / k^)Ai>  = T^i/Ai, 


(6.22) 


and  longitudinal  field, 


A — k^ki,/k^Au  — L^uAu 


(6.23) 


(on  which  the  action  does  not  depend),  satisfy 

dAT 

dL  - _t2AT  4-  nT 
- « X-n  + rj^ 


dr 


dA\i  L 

-d7  = r]» 


(6.24) 

(6.25) 


rj-i  T 

and  we  have  used  the  obvious  notation  r]j^  = 

Whereas  the  transverse  field  clearly  has  a Gaussian  equilibrium  distribu- 
tion, the  longitudinal  field  is  proportional  to  a pure  Brownian  motion,  and 
hence  has  no  equilibrium  distribution.  As  a consequence  of  this  fact  it  may  be 
shown  that  the  two  point  function 


E(ApAv) 


(6.26) 


does  not  converge  asr->  oo.  For  example,  if  we  take  as  an  initial  condition 


Afi(k,r  = 0)  = p<K^2) 

where  4>(k^)  is  a random  variable  with  a Gaussian  distribution 


(6.27) 


P{(f>)  oc  exp  J dn 


k^ 


(6.28) 
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and  average  over  this  distribution,  then  we  find  that  for  large  r,  the  propagator 
becomes 


6n(k  + q) 

(6.29) 

Whereas  the  term  proportional  to  r clearly  does  not  converge,  we  recognize 
the  finite,  r independent  part  as  the  a gauge  propagator.  Thus  the  choice 
of  initial  conditions  effectively  corresponds  to  a choice  of  gauge.  The  central 
point  is  that,  while  gauge  dependent  quantities  have  poorly  defined  equilibrium 
limits,  gauge  covariant  quantities,  such  as  E{FiiUFU)(T)  in  which  the  longitudinal 
field  does  not  appear,  do  converge.  In  this  sense  stochastic  quantization  does 
not  require  explicit  gauge  fixing.  We  remark,  however,  that  while  not  required, 
there  are  more  general  formulations  of  stochastic  quantization  which  do  fix  the 
gauge  by  considering  the  expanded  class  of  transformations 


E(A^k)Av{q))  ->  (2tt)" 


-j~2  — (1  — cx)kfiki//k'^  + ‘Irk^ki/ / k 


A^x,  t ) -»  An(x,  r)  + A(x,  r) 


(6.30) 


In  this  formalism,  the  longitudinal  field  will  also  have  an  equilibrium  distribu- 
tion, rendering  even  gauge  dependent  quantities  finite  in  the  limit  r — ► oo. 

Even  if  decidedly  formal  in  character,  the  functional  integral  approach  to 
stochastic  quantization  is  too  elegant  not  to  receive  at  least  a brief  description 
here.  Note  that  general  time  dependent  averages  are  given  by 

E(<t>(Ti)  . . . 4>{rn))  = J T>r/  </>(ti)  . . . (f>(Tn)  e~4  / dT1l  (r)  (6.31) 

and  hence  we  are  interested  in  the  properties  of  the  partition  function 


Z = J Vq  e~y  dTT]2^ 


(6.32) 
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Viewing  the  Langevin  equation  as  a change  of  variable  rj  — > <f>,  we  see  that  this 
gives  rise  to  the  Jacobian 


8t)(t) 


6<K?) 


d 82S\  c, 

i;  + sh)S(t-t) 


(6.33) 


which,  in  a manner  much  reminiscent  of  the  Faddeev-Popov  trick  for  inserting 
gauge  constraints,  may  be  formally  exponentiated  with  the  help  of  anticom- 
muting ghost  fields  c,  c 


jA=Jvcmexp(-J 


dr  c 


S2S 

dr  _r  82(j) 


+ 


(6.34) 


Inserting  unity  in  the  form 

'-h 


Sr/ 


8(f) 


r , d(t>  , 


(6.35) 


into  the  partition  function  and  integrating  out  the  Gaussian  variable  77,  we  see 
that 

(6.36) 


Z - J V(f>VcVc  e~S'ff  . 


Under  the  proper  assumptions,  the  resulting  effective  action 

2 / J 82S 


*.„  = /*  (s*g)  -*(s* 


82</> 


(6.37) 


can  be  shown  to  be  invariant  under  a supersymmetry  transformation 


6<j>  = ec  + ce 


8c  — e 


8c  — e 


8d4  _ 6S\ 

\dr  8(f)) 
}_d£  _ 8S\ 
\ dr  8(f) ) 


(6.38) 

(6.39) 

(6.40) 


The  significance  of  this  supersymmetry  has  been  much  discussed  in  the  litera- 
ture. For  example,  the  above  transformations  and  the  resulting  supersymmetry 
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are  essentially  the  inverse  of  the  Nicolai  map  [42]  for  supersymmetric  quan- 
tum mechanics.  Also  the  connection  with  the  dimensional  reduction  scheme 
proposed  by  Parisi  and  Sourlas  [43]  has  been  discussed.  There  are  formal 
arguments  indicating  that  the  supersymmetry  is  responsible  for  the  proper 
convergence  of  the  Langevin  process  to  the  desired  ground  state.  However,  as 
Klauder  and  Ezawa  [44]  have  pointed  out,  care  must  be  taken  in  imposing  the 
correct  boundary  conditions  on  all  fields.  In  demonstrating  the  supersymme- 
try one  partially  integrates  the  purely  bosonic  terms  in  the  effective  action  and 
drops  the  supposed  boundary  term 


/ 


dr 


d<t>8S 

dr  8(f) 


(6.41) 


However,  the  above  integral  should  in  fact  be  interpreted  as  a stochastic  inte- 
gral, and  hence  requires  a prescription  (for  example,  Ito  or  Stratonovich)  for 
its  definition. 

As  mentioned  in  the  beginning  of  this  chapter,  most  work  in  this  field  has 
presumed  a real  valued  action,  achieved  via  a Wick  rotation  to  Euclidean  space. 
However,  this  is  not  always  the  case,  and  studies  have  been  done  directly  in 
Minkowski  space.  Note  that  simulating  with  an  action  ihS  is  the  extreme 
example  of  complex  Langevin.  To  ensure  convergence  of  the  process  for  scalar 
fields  in  this  case,  a small  imaginary  mass  term  ^e<^>2  may  be  introduced  into 
the  action  S.  In  momentum  space  this  leads  to  the  free  equation 

^ = ilk?  — m2  + ie)<f>  + q (6.42) 

OT 

which  leads  to  the  two  point  function 


E(<f>(k,  r)(j)(q , r))  = (2?r )n8(k  + q)- ^ TT~ 

k + le 


(l-e2 i(k2~m2+ie)T^  . (6-43) 
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And  we  obtain  the  desired  Feynman  propagator  in  the  limit  r — > oo.  This 
procedure  in  Minkowski  space  has  also  been  studied  numerically  for  the  anhar- 
monic  oscillator  with  very  satisfactory  results  [16].  Complex  Langevin  may 
also  arise  in  applying  stochastic  quantization  to  gravity  and  to  general  tensor 
fields.  There,  a covariant  form  of  the  Langevin  equation  leads  to  non-positive- 
definite  noise  correlations,  which  may  be  formulated  in  terms  of  complex  valued 
fields. 

Finally  we  mention  an  equivalent  formulation  based  on  the  definition  of 
‘pseudo  white  noise’  [45].  In  this  scenario  the  Langevin  equation  comprises  of 
purely  real  drift  and  diffusion  terms.  However,  functional  averages  are  defined 
via  a complex  valued  distribution 

= JfVrl  exP  J dTrl2^j  (6-44) 

which  give  rise  to  the  complex  valued  expectation  values 

E(mvj)  = 2 iSij  . (6.45) 

This  formalism  has  been  explicitly  checked  to  give  the  usual  result  for  at  least 
quadratic  actions. 


CHAPTER  7 
CONCLUSIONS 


In  this  work  we  have  attempted  to  understand  and  further  develop  the 
underlying  theory  of  the  complex  Langevin  method  for  simulating  systems 
defined  by  a complex  valued  action.  This  method  has  an  innate  advantage 
over  standard  Monte  Carlo  (MC)  techniques  in  that  it  uses  the  entire  complex 
action  to  define  the  relevant  stochastic  process,  and  hence  a priori  avoids  any 
potential  sign  problem  which  often  results  in  MC  when  dividing  by  the  average 
of  the  complex  phase.  The  most  serious  limitation  of  the  CL  method  has 
been  its  occasional  failure  to  converge  to  the  desired  stationary  state,  either 
by  solutions  which  diverge  in  time,  or  by  long  time  averages  which  simply 
appear  to  converge  to  the  wrong  result.  The  proper  control  of  this  convergence 
in  practical  situations  remains  an  open  question,  even  though  the  method  of 
kernel  control  works  quite  well  for  some  model  problems. 

We  have  instead  addressed  the  more  modest  question  of  whether  it  is  possi- 
ble to  find  practical  numerical  criteria  which  will  at  least  determine  a posteriori 
whether  or  not  a given  simulation  has  indeed  converged  properly.  Toward  this 
end  we  have  investigated  the  general  conditions  under  which  a process  defined 
by  a complex  Langevin  equation  may  converge  to  the  proper  stationary  state. 
For  nonsingular  actions,  we  found  that  the  time  independence  of  the  ensemble 
expectation  values  is  necessary  and  sufficient.  This  was  demonstrated  first  on 
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compact  spaces  by  showing  that  expectation  values  have  unique  stationary  val- 
ues, in  spite  of  the  existence  of  ‘spurious  solutions’  to  the  pseudo  Fokker-Planck 
equation.  We  have  also  argued  this  for  the  noncompact  case. 

A further  useful  condition  is  found  from  the  Riemann-Lebesgue  lemma, 
which  requires  that  the  ensemble  expectation  values  of  the  characteristic-like 
function  E(e^z)  must  decrease  to  zero  as  |A:|  — ► oo.  Surprisingly,  we  find 
that  it  is  not  strictly  necessary  for  a pseudo  Fokker  Planck  equation  for  the 
complex  density  to  exist  in  the  sense  of  classical  function  theory;  indeed  it  may 
consistently  be  interpreted  as  an  equation  in  the  sense  of  distributions. 

In  the  model  problem  we  have  studied,  we  have  found  by  solving  a closed 
set  of  equations  for  the  moments  that  the  complex  Langevin  prescription  is 
actually  valid  much  more  often  than  numerical  integration  of  the  stochastic 
differential  equation  leads  one  to  believe.  We  conclude  that  this  is  due  to  the 
slow  convergence  of  the  process,  as  opposed  to  being  due  to  the  existence  of 
a degenerate  stationary  state.  We  also  found  that  the  mapping  from  a two 
dimensional  probability  density  to  a complex  valued  pseudo  density  is  far  from 
unique,  as  we  have  demonstrated  with  stochastically  similar  solutions.  For 
singular  actions,  we  are  yet  unable  to  prove  a general  result  on  the  convergence 
of  such  processes.  However,  we  did  find  empirically  on  model  problems  that 
ergodicity  is  a very  accurate  indicator  of  the  accuracy  of  the  simulation  results. 

The  theoretical  studies  described  above  were  greatly  motivated  by  the  de- 
sire to  apply  complex  Langevin  methods  to  the  problem  of  two-dimensional 
lattice  fermions.  For  many  such  interacting  systems,  there  are  presently  no 
reliable  numerical  studies  available,  as  standard  MC  techniques  are  hampered 
by  the  as  yet  unresolved  sign  problem.  Complex  Langevin  algorithms  for  inter- 
acting fermions  have  been  based  on  a spin  coherent  state  representation  of  the 
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path  integral.  Previously  complex  Langevin  simulations  were  effective  only  for 
one-dimensional  systems,  as  large  numerical  factors  resulting  from  the  asym- 
metric Jordan- Wigner  ordering  of  the  two-dimensional  lattice  were  present  in 
the  action.  This  problem  may  be  overcome  by  a suitable  choice  of  measure  on 
the  space  of  coherent  states,  and  we  have  found  this  numerically  to  be  quite 
successful.  While  the  results  of  complex  Langevin  for  such  systems  are  quite 
encouraging  for  the  free  fermion  case  we  studied  extensively,  the  process  is  not 
sufficiently  ergodic  to  reliably  compute  exact  results,  although  we  find  the  qual- 
itative behavior  to  be  well  reflected  in  the  numerical  simulations.  We  speculate 
that  this  lack  of  ergodicity  is  due  to  singular  nature  of  the  drift  terms  resulting 
from  the  vanishing  coherent  state  overlap  functions.  However,  it  is  encour- 
aging that  the  introduction  of  interaction  terms  to  the  fermion  Hamiltonian 
do  not  complicate  the  algorithm,  and  indeed  free  and  interacting  fermions  are 
treated  virtually  the  same  in  the  formalism.  Further  studies  and  refinements 
are  required  before  the  complex  Langevin  method  may  be  reliably  applied  to 
such  systems. 
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